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Introduction 

The  “Time-resolved  Spectral  Optical  Breast  Tomography”  research  project  aims  to  develop  a  near  real¬ 
time  three-dimensional  (3D)  spectral  tomographic  imaging  algorithm  with  use  of  the  cumulant 
approximation  to  radiative  transfer  to  model  light  propagation  in  tissues.  This  project  in  the  three  years' 
period  involves  theoretical  study  of  photon  migration  in  turbid  media,  implementing  an  enhanced  three- 
dimensional  (3D)  tomographic  imaging  algorithm,  and  performing  image  reconstruction  for  experimental 
data  obtained  from  tissue  phantoms  and  ex  vivo  tissues.  The  research  generated  14  journal  papers,  1  patent 
application,  and  23  conference  papers  and  presentations  during  this  period. 

Body 

The  tasks  performed  during  the  three  years'  period  include  theoretical  study  of  photon  migration  in  turbid 
media,  implementing  an  enhanced  3D  tomographic  imaging  algorithm,  and  performing  image 
reconstructions  for  experimental  data  obtained  from  tissue  phantoms  and  ex  vivo  tissues. 

Theoretical  study  of  photon  migration  in  turbid  media 

We  extended  the  photon  transport  model  for  light  migration  in  turbid  media  based  on  a  cumulant 
approximation  to  radiative  transfer  to  a  bounded  medium  with  planar  geometries[l]  (Task  1.1).  This 
extension  makes  the  cumulant  model  more  suitable  in  practical  applications  which  always  involve  finite 
geometries.  A  “reshaped  cumulant”  solution[2,  Appendix  1]  is  developed  to  improve  upon  the  previous 
second  order  cumulant  solution  to  radiative  transferfl]  described  by  a  Gaussian  distribution  in  two  aspects: 
1)  separating  the  ballistic  component  from  the  scattered  component  to  ensure  that  the  summation  in 
expressions  is  convergent;  and  2)  enforcing  the  causality  condition  to  ensure  that  no  light  travels  faster  than 
the  speed  of  light.  Time-resolved  profiles  obtained  using  the  analytical  form  were  compared  with  those 
obtained  by  the  Monte  Carlo  simulation,  for  both  transmission  and  backscattering.  A  significant 
improvement  of  the  accuracy  of  the  “reshaped  cumulant”  solution  over  the  original  cumulant  solution  is 
observed.  The  calculating  time  using  our  analytical  form  is  much  faster  than  that  using  the  Monte  Carlo 
approach,  usually  at  least  104  times  faster  (Task  1.1). 

We  also  developed  a  novel  Monte  Carlo  algorithm  (EMC)  [3,  Appendix  2]  to  simulate  polarized  light 
propagation  in  tissues  and  other  turbid  media  (Task  1.1).  This  Monte  Carlo  method  is  based  on  tracing  the 
multiply  scattered  electric  field  instead  of  the  conventional  Stokes  vector.  The  unique  feature  of  EMC 
makes  it  possible  to  simulate  also  coherent  properties  of  multiply  scattered  light.  The  algorithm  of  EMC  is 
straightforward  and  can  be  easily  adapted  to  simulate  the  propagation  of  polarized  light  in  optically  active 
or  gain  medium. 

Implementing  an  enhanced  3D  tomographic  imaging  algorithm 

We  enhanced  the  3D  tomographic  image  reconstruction  algorithms[4,  5]  by  using  the  new  cumulant 
transport  model  (Task  1.2)  and  making  use  of  spectral  information  when  observations  of  multiple 
wavelengths  are  available  (Task  1.3).  We  improved  the  3D  tomographic  image  reconstruction  algorithm  to 
use  a  L-curve  method  guided  by  the  signal-to-noise  ratio  of  the  dataset  to  determine  the  regularization 
parameter  (Task  1 .4). 

By  scanning  a  point  source  on  the  grids  of  the  input  plane  of  a  slab  and  measuring  light  intensities  on  a 
detector  array  on  the  exit  plane  of  the  slab,  a  set  of  four-dimensional  (4D)  data  is  formed.  The  spectral 
information  adds  an  additional  dimension  of  the  data.  The  optimal  approach  to  analyze  this  huge  dataset  is 
studied  (Task  1.5).  A  novel  optical  imaging  approach  based  on  independent  component  analysis  to  analyze 
such  massive  dataset  is  proposed  and  implemented[6-9;  Appendix  3,  Appendix  4], 
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The  image  reconstruction  algorithm  is  a  new  approach  for  optical  imaging  and  localization  of  objects  in 
turbid  media  that  makes  use  of  the  independent  component  analysis  (1CA)  from  information  theory. 
Experimental  arrangement  realizes  a  multisource  illumination  of  a  turbid  medium  with  embedded  objects 
and  a  multidetector  acquisition  of  transmitted  light  on  the  medium  boundary.  The  resulting  spatial  diversity 
and  multiple  angular  observations  provide  robust  data  for  three-dimensional  localization  and 
characterization  of  inhomogeneities  embedded  in  a  turbid  medium.  ICA  of  the  perturbations  in  the  spatial 
intensity  distribution  on  the  medium  boundary  sorts  out  the  embedded  objects,  and  their  locations  are 
obtained  from  Green's  function  analysis  based  on  any  appropriate  light  propagation  model. 

We  worked  with  the  experimental  group  at  the  Institute  of  Ultrafast  Lasers  and  Spectroscopy  and  to  test 
and  verify  our  image  reconstruction  algorithm  (Task  1.5,  Task  2.1,  Task  2.2).  We  have  successfully 
enhanced  the  algorithm  to  import  data  from  the  experimental  group,  to  locate  and  characterize 
inhomogeneities  within  the  turbid  medium,  and  to  generate  the  two  dimensional  (2D)  cross  sections  of  the 
inhomogeneities.  The  algorithm  is  tested  and  verified  to  work  well  for  absorption,  scattering,  and/or 
fluorescence  inhomogeneities  embedded  in  phantoms  using  different  wavelengths  (Task  2.1,  Task  2.2). 

To  access  the  efficacy  of  a  linear  inversion  scheme  in  image  reconstruction  of  human  breasts,  we  also 
studied  the  nonlinear  effect  of  the  multiple  passages  of  an  absorption  inhomogeneity  of  finite  size  deep 
inside  a  turbid  medium  on  optical  imaging  using  the  cumulant  solution  to  radiative  transfer  (Task  1).  We 
derived  the  analytical  nonlinear  correction  factor  which  agrees  excellently  with  the  predictions  from  the 
Monte  Carlo  simulations.  We  concluded  that  the  effect  of  the  nonlinear  multiple  passages  of  an  absorption 
site  on  optical  imaging  only  becomes  appreciable  when  the  size  of  the  inhomogeneity  reaches  1 0  times 
transport  mean  free  pat  or  larger  for  human  tissues[10,l  1]. 

Image  reconstructions  for  experimental  data 

We  have  successfully  performed  image  reconstructions  for  experimental  data  obtained  from  tissue 
phantoms  and  ex  vivo  tissues  with  absorption,  scattering,  and/or  fluorescence  inhomogeneities  embedded 
(Task  1 .5,  Task  2. 1 ,  Task  2.2,  Task  2.3).  To  optimize  the  algorithm  to  achieve  an  in  vivo  mammography, 
we  have  studied  the  performance  of  the  reconstruction  algorithm  on  breast  phantoms  and  ex  vivo  breast 
tissues  (Task  1.5,  Task  2.3).  A  total  of  about  10  cases  for  absorption,  scattering,  or  fluorescence 
inhomogeneities  embedded  in  phantoms  were  investigated.  In  each  case,  the  inhomogeneities  were  correctly 
located  and  characterized.  A  two-dimensional  (2D)  cross  section  map  of  the  inhomogeneities  was  then 
generated  afterwards  based  on  back-projection  [6-9,  Appendix  3,  Appendix  4].  One  notable  case  is  for  a 
breast  phantom  borrowed  from  University  College  London  which  has  a  thickness  of  55mm  (60  transport 
mean  free  path),  comparable  to  a  real  breast.  All  the  four  inhomogeneities  were  correctly  reconstructed  [6, 
Appendix  3],  including  the  one  of  the  lowest  contrast  which  has  only  10%  higher  scattering  than  the 
background  and  was  believed  to  be  undetectable  [12], 

Our  experimental  results  show  that  the  image  reconstruction  algorithm  is  superior  to  conventional  photon 
migration  reconstruction  algorithms  in  resolving  low  contrast  and  small  inhomogeneities  [6-9].  In  addition, 
this  approach  is  applicable  to  different  medium  geometries,  can  be  used  with  any  suitable  photon 
propagation  model,  and  is  amenable  to  near-real-time  imaging  applications.  This  approach  may  be  very 
useful  in  detecting  small  tumors  in  breasts  at  early  stage  of  development. 


Key  Research  Accomplishments 

•  Extended  cumulant  solution  of  radiative  transfer  to  planar  geometries  making  it  more  suitable 
for  practical  applications  which  involve  finite  boundaries. 

•  Developed  a  “reshaped  cumulant”  solution  for  photon  migration  in  tissues  and  other  turbid 
media,  improving  upon  the  previous  cumulant  solution  to  radiatve  transfer. 
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•  Developed  an  Electric  field  Monte  Carlo  (EMC)  algorithm  for  simulating  polarized  light 
propagation  in  turbid  media  whose  unique  feature  includes  simulating  coherence  phenomenon 
of  multiply  scattered  light. 

•  Developed  and  enhanced  the  3D  tomographic  image  reconstruction  algorithm  by  using  the 
cumulant  transport  model. 

•  Developed  the  criterion  of  the  optimal  regulation  parameter  for  inverse  image  reconstruction 
related  to  the  noise  presented  in  measurements  and  enhanced  the  3D  tomographic  image 
reconstruction  algorithm  to  use  a  L-curve  method  guided  by  the  signal-to-noise  ratio  of  the 
dataset  to  determine  the  regularization  parameter. 

•  Derived  the  nonlinear  correction  factor  of  multiple  passages  of  an  absorption  inhomogeneity  by 
a  photon  for  optical  imaging  and  provided  a  measure  of  the  efficacy  of  linear  inversion 
schemes  and  extended  the  range  of  applicability  of  the  linear  inversion  scheme  for  optical 
imaging. 

•  Developed  a  novel  3D  image  reconstruction  algorithm  “optical  imaging  using  independent 
component  analysis  (OPTICA)”. 

•  Tested  and  demonstrated  the  efficacy  of  OPTICA  in  imaging  absorption,  scattering,  and/or 
fluorescence  inhomogeneities  in  turbid  media,  in  particular,  for  low  contrast  small 
inhomogeneities. 


Reportable  Outcomes 
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Conclusions 

The  work  carried  out  during  the  three  years'  period  leads  to  the  following  conclusions.  First,  the  cumulant 
transport  model  provides  a  more  accurate  model  than  the  conventional  diffusion  model  for  the  description 
of  light  propagation  in  turbid  media  such  as  human  breasts.  Second,  the  “reshaped  cumulant”  solution 
improves  over  the  previous  cumulant  solution,  providing  a  even  more  accurate  model  for  light  propagation 
in  turbid  media.  Third,  the  optimal  regularization  of  image  reconstruction  depends  on  the  noise  presented  in 
the  measurements;  proper  modeling  of  the  noise  and  appropriate  regularization  improves  the  quality  of 
image  reconstruction.  Fourth,  the  correction  for  the  nonlinear  effect  of  the  multiple  passages  of  an 
absorption  site  by  light  was  shown  to  be  essential  in  optical  imaging  to  characterize  properly 
inhomogeneities  strong  in  absorption.  Fifth,  a  novel  image  reconstruction  algorithm  “optical  imaging  using 
independent  component  analysis  (OPTICA)”  has  been  developed,  able  to  resolve  low  contrast  small 
absorption,  scattering,  and/or  fluorescence  inhomogeneities  in  turbid  media.  Sixth,  the  theoretical  formalism 
and  computer  algorithm  for  3D  tomographic  image  reconstruction  shows  with  experimental  data  the 
potential  to  provide  fast  3D  images  of  small  tumors  in  breasts  at  early  stage  of  development. 
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An  analytical  expression  of  the  particle  distribution  based  on  an  analytical  cumulant  solution  of  the  time- 
dependent  elastic  Boltzmann  transport  equation  (BTE)  is  presented.  This  expression  improves  upon  the  previ¬ 
ous  second  order  cumulant  solution  of  the  BTE  described  by  a  Gaussian  distribution  in  two  aspects:  (1) 
separating  the  ballistic  component  from  the  scattered  component  to  ensure  that  the  summation  in  expressions 
is  convergent;  and  (2)  enforcing  the  causality  condition  to  ensure  that  no  particle  travels  faster  than  the  free 
speed  of  the  particles.  Time-resolved  profiles  obtained  using  the  analytical  fonn  are  compared  with  those 
obtained  by  the  Monte  Carlo  simulation,  for  both  transmission  and  backscattering.  The  calculating  time  using 
our  analytical  form  is  much  faster  than  that  using  the  Monte  Carlo  approach. 

DOT:  10.1 103/PhysRevE.7 1.04 1202  PACS  number(s):  42.25.Dd,  05.60.Cd,  05.20.-y,  42.25.Fx 


I.  INTRODUCTION 

The  time-dependent  elastic  Boltzmann  transport  equation 
(BTE)  describes  the  particle  (and  light,  acoustic  wave,  etc.) 
propagation  with  time  in  a  scattering  medium,  where  the 
particles  suffer  multiple  scattering  by  randomly  distributed 
scatterers.  The  BTE  is  also  called  the  radiative  transfer  equa¬ 
tion  in  light  propagation  [1-3].  The  solutions  of  the  elastic 
BTE  are  applied  in  broad  areas,  such  as  atmospheric  science, 
medical  imaging,  and  solid  state  physics. 

An  example  is  the  approach  to  optical  imaging  of  human 
tissue  that  is  often  called  “diffusion  tomography,”  because 
the  theoretical  model  is  built  based  on  the  solution  of  the 
diffusion  equation.  The  diffusion  equation  is  the  lowest  order 
approximation  of  the  radiative  transfer  equation,  which  has 
significant  error  when  the  distance  between  a  voxel  and  a 
source  is  short.  However,  the  contribution  from  these  voxels 
near  the  source  to  the  measured  signals  is  much  larger  than 
that  from  voxels  deep  inside  body.  Hence,  for  accurate  im¬ 
aging  the  theoretical  model  should  be  based  on  solution  of 
the  radiative  transfer  equation.  A  similar  procedure  can  be 
applied  to  images  of  cloud  distribution  obtained  using  a  lidar 
arranged  on  a  satellite,  which  requires  knowledge  of  the  mul¬ 
tiple  scattering  effect  of  water  drops  distributed  in  the  cloud 
on  the  time-resolved  backscattering  signals.  In  both  ex¬ 
amples,  the  size  of  the  scatterers  can  be  nearly  equal  to  or 
larger  than  the  wavelength  of  light,  leading  to  a  large  aniso¬ 
tropic  factor.  The  use  of  low-frequency  sound  to  detect  oil¬ 
bearing  layers  deep  under  the  ocean  floor  is  another  example. 

Currently,  numerical  approaches,  including  Monte  Carlo 
simulations,  are  the  main  methods  in  solving  the  BTE  [4-6], 
Numerical  solution  of  the  BTE  is  a  cumbersome  task,  since 
the  particle  distribution  7(r,s,/)  is  a  function  of  position  r, 
angle  s,  and  time  /,  in  a  six-dimensional  space  of  parameters. 
An  analytical  expression  for  7(r,s,/)  with  quantitative  accu¬ 
racy  can  greatly  reduce  the  computation  burden  in  modeling 
particle  and  light  propagation  in  scattering  media,  which  is 
essential  for  imaging  in  turbid  media,  because  the  inverse 


reconstruction  process  calls  the  forward  model  many  times. 

Recently,  we  have  developed  an  analytical  solution  of  the 
time-dependent  elastic  BTE  in  an  infinite  uniform  medium 
with  an  arbitrary  phase  function  [7,8],  The  exact  spatial  cu- 
mulants  of  7(r,s,/)  up  to  an  arbitrary  high  order  at  any  angle 
and  any  time  have  been  derived.  A  cutoff  at  second  order  of 
the  cumulants  7(r,s,/)  can  be  approximately  expressed  by  a 
Gaussian  distribution,  which  has  the  exact  first  cumulant  (the 
position  of  the  center  of  the  distribution)  and  the  exact  sec¬ 
ond  cumulant  (the  half-width  of  the  spread  of  the  distribu¬ 
tion).  The  cumulant  solution  of  BTE  has  been  extended  to 
the  case  of  a  polarized  photon  distribution,  and  to  semi¬ 
infinite  and  slab  geometries.  Using  a  perturbation  method, 
the  distribution  7(r,s,/)  in  a  weak  heterogeneous  medium 
can  be  calculated  based  on  the  cumulant  solution  of  the  BTE. 

The  analytical  cumulant  solution  of  the  BTE  obtained, 
although  it  has  exact  center  and  half-width,  is  not  satisfac¬ 
tory  in  two  respects.  First,  one  cannot  ensure  that  the  sum¬ 
mation  over  /  in  the  expressions  shown  in  Sec.  II  is  conver¬ 
gent  at  very  early  times.  Second,  a  remarkable  fault  of  the 
Gaussian  distribution  at  early  times  is  that  particles  at  the 
front  edge  of  the  distribution  travel  faster  than  the  free  speed 
of  the  particles  in  the  medium,  thus  violating  causality,  espe¬ 
cially  for  those  particles  moving  along  near  forward  direc¬ 
tions.  The  Gaussian  distribution  is  accurate  at  long  times  and 
in  the  backscattering  case,  since  many  collisions  lead  to  a 
Gaussian  distribution  according  to  the  central  limit  theorem. 

In  this  paper,  the  analytical  cumulant  solution  of  the  BTE 
has  been  improved  compared  to  our  previous  work  [7]  in 
these  two  respects.  For  solving  the  first  problem,  we  make  a 
separation  of  the  ballistic  component  from  the  total  7(r,s,/) 
and  compute  the  cumulants  for  the  scattered  component 
x(r,s,t)-  This  treatment  ensures  convergent  summation  over 
/.  Also  this  separation  provides  a  clearer  picture  of  particle 
propagation.  In  the  time-resolved  transmission  profile  the 
ballistic  component  is  described  by  a  sharp  jump  exactly  at 
the  ballistic  time,  separated  from  the  later  scattered  compo¬ 
nent.  For  solving  the  second  problem  two  approaches  are 
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used.  The  first  method  is  to  calculate  the  distribution  includ¬ 
ing  the  higher-order  cumulants,  based  on  our  work  in  Ref. 
[8],  However,  computation  of  high-order  cumulants  is  a 
cumbersome  task.  In  the  second  method  the  Gaussian  distri¬ 
bution  is  replaced  by  a  different-shaped  form,  which  satisfies 
causality,  and  maintains  the  correct  center  position  and  the 
correct  half-width  of  the  distribution  computed  by  our  ana¬ 
lytical  approach.  There  are  infinite  choices  of  the  shapes  of 
the  distribution  satisfying  these  conditions;  we  choose  a 
simple  analytical  form.  At  long  times,  the  reshaped  distribu¬ 
tion  tends  to  the  Gaussian  distribution.  Our  results  show  that 
the  reshaped  distribution  matches  that  obtained  using  Monte 
Carlo  simulation  much  better  than  the  Gaussian  distribution. 

The  paper  is  organized  as  follows.  In  Sec.  II  we  briefly 
review  the  main  results  of  the  analytical  cumulant  solution  of 
the  BTE.  Section  III  presents  a  separation  of  the  ballistic 
component  from  the  scattered  component,  which  makes  the 
summation  over  /  convergent.  Section  IV  improves  the  dis¬ 
tribution  at  early  times  using  two  approaches,  and  presents 
the  numerical  result  compared  with  the  Monte  Carlo  simula¬ 
tion.  Section  V  is  devoted  to  discussion  and  conclusions. 


II.  THE  ANALYTICAL  CUMULANT  SOLUTION 
OF  THE  BTE 

The  elastic  Boltzmann  kinetic  equation  of  particles,  with 
magnitude  of  velocity  v,  for  the  distribution  function  7(r,s,f) 
as  a  function  of  time  t,  position  r,  and  direction  s,  in  an 
infinite  uniform  medium,  from  a  point  pulse  light  source, 
<S(r-r0)<5{s-So)<S(!-/o)>  is  given  by 

dl(r,s,t)/dt  +  vs  ■  Vr/(r,s,t)  +  /xj(r,s,t) 

M'Wr, »',-)*• -Mr.M) 

+  <5(r-r0)<5(s-s0)<S(r-ro)  (1) 

where  jxs  is  the  scattering  rate,  /xa  is  the  absorption  rate,  and 
P( s,s')  is  the  phase  function,  normalized  to  //>(s,s')cfs'=  1. 
The  phase  function  is  assumed  to  depend  only  on  the  scat¬ 
tering  angle  in  an  isotropic  medium.  Under  this  assumption, 
an  arbitrary  phase  function  can  be  handled.  We  expand  the 
phase  function  in  Legendre  polynomials  with  constant  coef¬ 
ficients, 


P(s,s')=y-^Ja,Pl{s-s').  (2) 

47T  i 

Recently,  we  have  developed  a  different  approach  to  ob¬ 
tain  an  analytical  solution  of  the  BTE  in  an  infinite  uniform 
medium,  based  on  a  cumulant  expansion  [7,8]. 

We  briefly  review  the  concept  of  the  “cumulant”  in  a  one¬ 
dimensional  (ID)  case.  Consider  a  random  variable  x,  with  a 
probability  distribution  function  /(x).  Instead  of  using /(x)  to 
describe  the  distribution,  we  define  the  «th  moment  of  x, 

(3) 

and  correspondingly  the  «th  cumulant  (x")c  defined  by 


(00  \  cc 

2  (x”)c(;7)"/n! )  =  (exp(frx))  =  2  (x")(it)"/n\  (4) 

The  first  cumulant  (x)c  is  the  mean  position  of  x.  The 
second  cumulant  (x2)c  represents  the  half-width  of  the  distri¬ 
bution.  The  higher  cumulants  are  related  to  the  detailed 
shape  of  the  distribution.  For  example,  (x3)c  describes  the 
skewness  or  asymmetry  of  the  distribution,  and  (x4)c  de¬ 
scribes  the  “kurtosis”  of  the  distribution,  that  is,  the  extent  to 
which  it  differs  from  the  standard  bell  shape  associated  with 
the  Gaussian  distribution  function.  The  cumulants  hence  de¬ 
scribe  the  distribution  in  an  intrinsic  way  by  subtracting  off 
the  effects  of  all  lower-order  moments.  In  the  3D  case,  the 
first  cumulant  has  three  components,  the  second  cumulant 
has  six  components,  and  so  on. 

We  have  derived  an  explicit  algebraic  expression  for  the 
spatial  cumulants  at  any  angle  and  any  time,  exact  up  to  an 
arbitrarily  high  order  n  [8],  This  means  the  distribution  func¬ 
tion  l(r,s,t)  can  be  computed  to  any  desired  accuracy.  At  the 
second  order,  n= 2,  an  analytic  explicit  expression  for  distri¬ 
bution  function  7(r,s,/)  is  obtained  [7,8].  This  distribution  is 
Gaussian  in  position,  which  is  accurate  at  later  times,  but 
only  provides  the  exact  mean  position  and  the  exact  half¬ 
width  at  early  times. 

The  Gaussian  distribution  of  the  second-order  cumulant 
solution  is  written  as 


7(r,s  ,t) 


F{$,SqJ)  1 
(4ir)3/2  (detS),/2 


Xexp 


~{B  ])ap(r-rc)a(r-rc)p 


(5) 


where  F(s,s0,l)  is  the  total  angular  distribution  F(s,s0,t) 
=fl(r,s,t)dr,  which  has  the  following  exact  expression: 


„  21+  1 

F(s,s0,t)  =  exp(-  iiJ)2j  — —  exp(-g/t)P/(s  •  s0) 

,  4  7T 

=  exp(-  /V)2  exp(“g//)S  YiJs)Y*/m(s0), 


where 

=  “«/( 2/+1)].  (7) 

Two  special  values  of  gt  are  g0=0,  which  follows  from  the 
normalization  of  7>(s,s'),  and  g]=u//tr,  where  7tr  is  the  trans¬ 
port  mean  free  path,  defined  by  /tr=u/[/u(l-(cos  0))], 
where  (cos  9)  is  the  average  of  ss'  with  P(s,s')  as  weight. 
In  Eq.  (6),  T/„,(s)  are  spherical  harmonics  normalized  to 
4i7/(2/+l). 

The  center  of  the  packet  (the  first  cumulant),  denoted  by 
rc,  is  located  at 


r\  =  G2  A,P,{cos  9)[(l  +  1  )f(g, - g/+I)  +  //(g, - g/-i)], 
/ 


(8) 
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rcx=G 2  A,P(/l)( cos  0)cos  <p[fig,- gh\) ~Agi~gM )L 

i 

(9) 


where  G=v  exp(-fiat)/F{s,s0,t),  ^/=(1 /4Tr)exp(-g/?),  and 
for  any  variable  x, 

/(.t)  =  [exp(xt)  -  1  ]/x .  (10) 


rc  is  obtained  by  replacing  cos  cp  in  Eq.  (9)  by  sin  <p.  In  Eqs. 
(8)  and  (9),  P\"'\ cos  6)  is  the  associated  Legendre  Emotion. 

The  square  of  the  average  spread  half-width  (the  second 
cumulant)  is  given  by 


Bap-vGAap- 


■rcjy  2, 


(ID 


where  all  the  coefficients  are  Einctions  of  angle  and  time: 


A„  =  X  A,P,{ cos  8) 


/(/-1)4ni  (/+1)(/+2)42) 


21-  1 


e\2)  + 


2/-1 

(/+  l)2 


21+3 


21+3 


E\4) 


(12) 


=  S  -A,P,{ COS  0) 


/(/-l)r,„  (/+!)(/+  2)  (2) 


2/-1 


21  +  3 


Kl-  D^m.  (/+D(/  +  2), 


+  - - -£T'  + 

21- 1  ' 


,(4) 


2/+  3 


±  2  ~ A \cqs  0)cos{2(p) 
;  2 


X 


1  1  ’  +  TT^-£l2)  -  E\ 


O)  _ . 
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21-  1 


2/  +  3 


21- 1  1  21+3 


e\4) 


(13) 


where  (+)  corresponds  to  Axx  and  (-)  corresponds  to  Ayy, 

Axy=A yx  =  2  ~A,P(2)(cos  8) sin(2<£) 

/  z 


1 


rr(2)  _  . 


1 


21+3  1  21-  1 


r(3). 
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21+3 


21-  1 


7(4) 


(14) 


FIG.  I .  The  moving  center  Rcz  and  the  diffusion  coefficients  Dzz 
and  Dxx  of  the  particle  density  function  as  functions  of  time  t. 


The  second  order  cumulant  approximation  for  the  particle 
density  distribution  N(r,t)  has  a  Gaussian  shape: 

(z-Rcf 
exp 


1 


1 


"(r’')  =  (4,-Z> 

::ut),/2  4ttD 

(x2+y2) 

Xexp 

4  Dxxvt 

with  a  moving  center  located  at 

AD„vt 


exp(-/v), 


7?J  =  t;[l  -exp(-g,?)]/g. 


(20) 


(21) 


and  the  corresponding  time-dependent  diffusion  coefficients 
are  given  by 


v  1  t 

D=-\ - 


3gl~g2 


[1  -exp(-g,t)] 


3Hgl  gl(gl“g2) 

2  3 

+  - :[1  -  exp(-g2t)]  -  — j[l  -  exp(-g|/)]2 

g2(g]-g2)  2  g, 


(22) 


A*r- A;.t-  X  2rt/P$°(cos  6>)cos  </j  2/  j 

2(/  +  2V+  1  E^+  1  £<4) 

21+3  1  2/-1  1  21+3  1  J 

(15) 

Ayz  is  obtained  by  replacing  cos  <p  in  Eq.  (15)  by  sin  (p.  In 
Eqs.  ( 1 2)— (1 5)  £'|i_4)  are  given  by 

£/')  =  [f(gi-g!-i)  -/(g/~g/-i)]/(g/-i -gi-i). 

(16) 

E{P  =  [figl  -  g/+2)  -f(g!  -  g/+ 1  )]/(g/+ 1  -  g/+2) , 

(17) 

£/3)  =  [figi — g/-i)  —  f\Kgi  -  g/- 1 ) , 

(18) 

e\  4)  =  [fig/  -  g/+ 1 )  -  /]%/  -  g/+ 1)  • 

(19) 

D=DV  =—  |  “ — I - r — — - 

“  3/  [g,  g2(g,-g2) 


[1  -exp(-g,/)] 


1 


g2(g\~g2) 


[1  -exp(-g2/)]  f. 


(23) 


Each  distribution  in  Eqs.  (5)  and  (20)  describes  a  particle 
“cloud”  anisotropically  spreading  from  a  moving  center,  with 
time-dependent  difEision  coefficients.  As  shown  in  Fig.  1,  at 
early  time  t—* 0,  the  mean  position  of  the  photons  moves 
along  the  s0  direction  with  speed  v,  and  the  diffusion  coeffi¬ 
cient  tends  to  zero.  These  results  present  a  clear  picture  of 
nearly  ballistic  motion  at  t—>0.  With  increase  of  time,  the 
motion  of  the  center  slows  down,  and  the  diffusion  coeffi¬ 
cients  increase  from  zero.  This  stage  of  particle  migration  is 
oEen  called  a  “snakelike  mode.”  At  late  times,  the  total  an¬ 
gular  distribution  function  tends  to  become  isotropic.  The 
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1(PS) 


FIG.  2.  Light  distribution  in  an  infinite  uniform  medium  as  a 
function  of  time  at  different  received  angles,  using  the  second  cu- 
mulant  solution  of  the  BTE  (Gaussian  distribution),  where  the  de¬ 
tector  is  located  at  z=10  mm  from  the  source  in  the  incident  direc¬ 
tion.  The  parameters  for  this  calculation  are  /„= 2  mm,  the 
absorption  length  /„=300  mm,  the  phase  function  is  computed  us¬ 
ing  Mie  theory  for  polystyrene  spheres  with  diameter  d=  1.11  fim 
in  water,  and  the  wavelength  of  the  laser  source  \  =  625  nm,  which 
gives  the  g  factor  g= 0.926. 


/A)(r,s,f)  =  exp[-  (fis  +  fia)t\S(r  -  vts0)S(s  -  s0).  (24) 

The  moments  of  the  ballistic  component  can  be  easily  calcu¬ 
lated.  When  s0  is  along  z,  we  have 

j  z”/b\r,s,t)d3r=  exp[(-  fis  +  /na)t](vt)"S(s  -  s0), 

(25) 

and  other  moments  related  to  z"1x"2/'3  (w2,«3  +  0)  are  zero. 

The  total  distribution  is  the  summation  of  the  ballistic 
component  and  the  scattered  component: 

7(r,s,/)  =  /A)(r,s ,/)  +  7(s)(r,M);  (26) 

hence,  the  moments  of  the  scattered  component  can  be  ob¬ 
tained  by  subtracting  the  corresponding  ballistic  moments 
from  the  moments  of  7(r,s,t).  For  example,  we  have 

J  zn/s\r,s,()d3r  =  J  znl(r,s,t)d3r-  J  z"7(4)(r,s ,i)d3r. 

(27) 

Notice  that 

<5(s-*o)  =  2[(2/+  I)/4t7]P/(s  •  s0).  (28) 


particle  density,  at  t>ltr/v  and  r>!u,  tends  toward  the 
center-moved  (l/tr)  diffusion  solution  with  the  diffusive  co¬ 
efficient  IJ 3.  Therefore,  our  solution  quantitatively  de¬ 
scribes  how  particles  migrate  from  nearly  ballistic  motion,  to 
snakelike  motion,  and  then  to  diffusive  motion. 

Figure  2  shows  the  calculated  distribution  as  a  function  of 
time  at  different  receiving  angles  in  an  infinite  uniform  me¬ 
dium,  computed  by  the  second  order  cumulant  solution, 
where  the  detector  is  located  at  5/tr  from  the  source  in  the 
incident  direction  of  the  source.  Figure  2  shows  anisotropy 
of  the  distribution  at  a  distance  of  5/,r  from  the  source.  This 
type  of  distribution  has  been  demonstrated  by  time-resolved 
experiments  [9]. 

The  analytical  solution  obtained,  although  it  has  the  exact 
center  and  half-width,  is  not  satisfactory  in  two  respects. 
First,  at  very  early  times,  exp(-g//)— > 1  for  all  /;  hence,  one 
cannot  ensure  that  summation  over  /  is  convergent.  Second, 
particles  at  the  front  edge  of  the  Gaussian  distribution  travel 
faster  than  the  speed  v,  thus  violating  causality. 


III.  SEPARATING  THE  BALLISTIC  COMPONENT  FROM 
THE  SCATTERED  COMPONENT 

In  order  to  make  the  summation  over  /  convergent,  we 
separate  the  ballistic  component  from  the  total  7(r,s,t),  and 
compute  the  cumulants  for  the  scattered  component  7*~s) 
X(r,s,/). 

The  ballistic  component  is  the  solution  of  the  homoge¬ 
neous  Boltzmann  transport  equation,  which  is  the  transport 
equation,  Eq.  (1),  without  the  “scattering  in”  term  [the  first 
term  in  the  right  side  of  Eq.  (1)].  The  solution  of  the  ballistic 
component  is  given  by 


Substituting  Eqs.  (25)  and  (28)  into  Eq.  (27),  the  correspond¬ 
ing  cumulants  for  the  scattered  component  fo\r,s,t)  can  be 
easily  obtained;  they  replace  Eqs.  (6),  (8),  and  (12)  by 

21+  1 

7^j)(s,s0>?)  =  exp(-  fiat)Z  1 — [exp(-  g,t) 
l  47 7 

-exp(-/y)]P/(s-s0),  (29) 


rfs)  =  G2  7j/(cos  8)-~-{exp{- g,t)[{l  +  l)/(g/-g/+i) 

/  477 

+  1)]  “  (2/  +  l)t  exp(-  /y)},  (30) 


=  2  P /(cos  e)~\  exp  (-git) 


(/+  D(/  +  2) 


21+3 
t\2l+  1) 


£;2)  +  ■ 


21-\E‘ 

(/+  l)2, 


21-  1 


£<3)  + 


■(4) 


21  +  3 


exp(—  /xst)  f . 


(31) 


The  expressions  for  the  other  components  of  the  first  and 
second  cumulants  are  unchanged,  provided  all  F{s,s0,t)  in  G 
in  Sec.  II  are  replaced  by  F^s\s,s0,t).  Note  that  Eq.  (28) 
actually  is  equal  to  zero  at  s^Sq,  and  there  is  no  ballistic 
component  in  these  directions. 

The  replacement  of  equations  in  Sec.  II  by  Eqs.  (29)— (3 1) 
greatly  improves  the  calculation  of  cumulants  at  very  early 
times.  By  the  subtraction  introduced  above,  the  terms  for 
large  /  approach  zero,  and  summation  over  /  becomes  con¬ 
vergent  at  very  early  times.  When  f— » 0,  g/=/x4.[l -a//(2/ 
+  1)]  [see  Eq.  (7)]  approaches  fis  for  large  /,  f[grgi±\)~t 
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FIG.  3.  Time-resolved  profile  of  the  backscattered  (180°)  photon 
intensity  inside  a  disk  with  center  at  r=0,  radius  R=\lu,  thickness 
dz=0.1/tr,  and  the  received  solid  angle  d cos  0=0.001,  normalized 
to  inject  one  photon.  The  Heyney-Greenstein  phase  function  with 
g=0.9  is  used,  and  1  /  /„ = 0 .  The  solid  curve  is  for  the  second  cu- 
mulant  solution  (Gaussian  distribution),  and  dots  are  for  the  Monte 
Carlo  simulation.  The  inset  diagram  shows  the  same  result  drawn 
using  a  logarithmic  scale  for  intensity. 

[see  Eq.  (10)],  and  7s{H4)  — 12/2  [see  Eqs.  (16)— (19)],  which 
results  in  cancellation  in  the  summand  for  large  /  at  very 
early  times. 

An  example  of  successful  use  of  this  replacement  is  the 
calculation  of  backscattering.  When  0=180°,  P/( cos  0)=1  or 
-1,  depending  on  whether  /  is  even  or  odd.  The  computed  rcz 
at  very  early  times  using  Eq.  (8)  oscillates  with  a  cutoff  of  I. 
But  the  computed  at  very  early  times  using  Eq.  (30) 
becomes  stable.  Calculation  shows  that  r^<1,=0  at  any  time 
when  0=180°. 

Figure  3  shows  the  computed  time  profile  of  the  back- 
scattering  intensity  7w(r,s ,/)  at  a  detector  centering  at  r=0 
and  received  at  an  angle  0=  1 80°,  which  actually  is  the  total 
backscattering  intensity  I(r,s,t)  because  s  +  s0,  compared 
with  the  Monte  Carlo  simulation.  The  absolute  value  of  the 
intensity,  as  well  as  the  shape  of  the  time-resolved  profile, 
computed  using  our  analytical  cumulant  solution  of  the  BTE 
match  well  with  those  of  the  Monte  Carlo  simulation.  The 
inset  diagram  is  the  same  result  drawn  using  a  logarithmic 
scale  for  intensity.  Note  that  this  result  of  backscattering, 
based  on  the  solution  of  the  BTE,  is  for  a  detector  located 
near  the  source,  different  from  other  backscattering  results 
based  on  the  diffusion  model,  which  is  only  valid  when  de¬ 
tector  is  located  at  a  distance  of  several  lu  from  the  source. 

IV.  SHAPE  OF  THE  PARTICLE  DISTRIBUTION 

If  the  cumulants  with  order  n  >  2  are  assumed  all  zero,  the 
distribution  becomes  Gaussian.  The  Gaussian  distribution  is 
accurate  at  long  times.  At  early  times,  particles  at  the  front 
edge  of  the  distribution  travel  faster  than  the  free  speed  of 
the  particles,  thus  violating  causality,  especially  for  particles 


FIG.  4.  Time-resolved  profile  of  transmitted  light  in  an  infinite 
uniform  medium,  computed  using  the  tenth  order  cumulant  solution 
(solid  curve),  the  second  order  cumulant  solution  (dotted  curve), 
and  the  diffusion  approximation  (thick  dots  curve),  compared  with 
that  of  the  Monte  Carlo  simulation  (discrete  dots).  The  detector  is 
located  at  z=6/lr  from  the  source  along  the  incident  direction,  and 
the  received  direction  is  0=0.  The  Heyney-Greenstein  phase  func¬ 
tion  with  g=0.9  is  used,  and  the  absorption  coefficient  1  //„=0. 

that  move  along  near  forward  directions.  In  the  following, 
two  approaches  are  used  for  overcoming  this  fault:  (A)  in¬ 
cluding  higher  cumulants;  and  (B)  introducing  a  reshaped 
distribution. 

A.  Calculation  including  high-order  cumulants 

We  have  performed  calculations  including  the  higher- 
order  cumulants  to  obtain  a  more  accurate  shape  of  the  dis¬ 
tribution.  The  Codes  for  calculation  are  designed  based  on  a 
formula  derived  in  Ref.  [8], 

Figure  4  shows  7(r,s,/)  with  the  detector  located  at  z 
=  6/tr  in  front  of  the  source  and  receiving  direction  along  0 
=  0,  computed  using  the  analytical  cumulant  solution  up  to 
tenth  order  of  the  cumulants  (solid  curve),  up  to  the  second 
order  cumulants  (dotted  curve),  in  the  diffusion  approxima¬ 
tion  (thick  dotted  curve),  and  the  Monte  Carlo  simulation 
(discrete  dots).  The  figure  shows  that  the  tenth  order  cumu¬ 
lant  solution  is  located  in  the  middle  of  the  data  obtained  by 
the  Monte  Carlo  simulation,  and  7(r,s,?)~0  before  the  bal¬ 
listic  time  tb=6ln/v.  The  second  order  cumulant  solution  has 
nonzero  7(r,s,t)  before  th,  which  violates  causality.  The 
computed  N{r,t)/4Tr  using  the  diffusion  model  has  a  large 
discrepancy  with  the  Monte  Carlo  simulation,  and  the  diffu¬ 
sion  solution  has  more  nonzero  components  before  tb,  which 
violates  causality. 

Using  the  second  order  cumulant  solution,  the  distribution 
function  can  be  computed  very  fast.  The  associated  Legendre 
functions  can  be  quickly  computed  using  recurrence  relations 
with  accuracy  limited  by  the  computer  machine  error.  It 
takes  1  min  to  produce  105  data  points  of  7(r,s,/)  on  a  per¬ 
sonal  computer.  On  the  other  hand,  in  order  to  reduce  the 
statistical  fluctuation  to  the  level  shown  in  Fig.  4,  109  events 
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are  counted  in  the  Monte  Carlo  simulation,  which  takes  tens 
of  hours  computation  time  on  a  personal  computer.  Compu¬ 
tation  of  high-order  cumulants  also  is  a  cumbersome  task, 
because  the  number  of  involved  terms  rapidly  grows  with 
increase  of  the  order  n.  Also,  for  a  distribution  function  that 
is  positive  definite,  the  Bust  theorem  states  that  the  existence 
of  nonzero  cumulants  of  any  order  higher  than  2  will  be 
accompanied  by  nonzero  cumulants  of  all  orders  [10].  There¬ 
fore,  no  matter  how  the  cutoff  at  a  finite  high  order  n  is 
taken,  the  cumulant  solution  of  the  BTE  cannot  be  regarded 
as  exact. 

B.  Reshaping  the  particle  distribution 

For  practical  applications,  we  use  a  semiphenomenologi- 
cal  model.  The  Gaussian  distribution  is  replaced  by  a 
different-shaped  form,  which  maintains  the  correct  center  po¬ 
sition  and  the  correct  half-width  of  the  distribution.  This  dis¬ 
tribution  satisfies  causality,  namely,  7(r,s,/)=0  outside  the 
ballistic  limit  vt.  There  are  an  infinite  number  of  choices  of 
shape  of  the  distribution  under  the  above  conditions.  We 
choose  a  simple  analytical  form  as  discussed  later.  At  long 
times,  the  half-width  of  the  distribution  a~(4B)U2,  with  B 
shown  in  Eq.  (11),  spreads  as  tm\  hence,  cr<vt  at  large  t, 
and  the  Gaussian  distribution  at  long  times  with  half-width  cr 
can  be  regarded  as  completely  inside  the  ballistic  sphere.  The 
reshaped  distribution  of  7(r,s,/)  hence  should  approach  the 
Gaussian  distribution  at  long  times. 

1.  ID  density 

We  first  discuss  the  one-dimensional  density  as  an  ex¬ 
ample  to  explain  our  model. 

The  Gaussian  distribution  of  ID  density  is  described  by 

m  =  (4-7r£).-i;?)~1/2exp[-  (z  -  7?f)2/(4ZL.o/)],  (32) 

where  Rcz  and  D..  are  given  in  Eqs.  (21)  and  (22).  As  shown 
in  Fig.  5,  although  the  ID  Gaussian  spatial  distribution  (the 
dashed  curve)  at  time  t=21xr/v,  Eq.  (32),  has  the  correct  cen¬ 
ter  and  half-width,  the  curve  deviates  from  the  distribution 
computed  by  the  Monte  Carlo  simulation  (dots),  and  a  re¬ 
markable  part  of  the  distribution  appears  outside  the  ballistic 
limit  vt=2lXr  At  early  times  the  spatial  distribution  is  not 
symmetric  about  the  center  Rc.  When  Rc  moves  from  the 
source  toward  the  forward  side,  causality  prohibits  particles 
appearing  beyond  vt.  This  requires  the  particles  in  the  for¬ 
ward  side  to  be  squeezed  in  a  narrow  region  between  Rc  and 
vt.  For  a  balance  of  the  parts  of  the  distribution  in  the  for¬ 
ward  and  backward  sides  of  Rc,  the  peak  of  the  distribution 
should  move  to  a  point  at  the  forward  side  and  the  height  of 
the  peak  should  increase.  The  earlier  the  time  t,  the  closer  is 
Rc  to  vt,  and  the  asymmetry  of  the  distribution  becomes 
stronger.  Based  on  this  observation  we  propose  the  following 
analytical  expression:  (1)  to  move  the  peak  position  of  the 
distribution  from  Rcz  to  zc,  where  the  parameter  zc  will  be 
determined  later;  (2)  to  take  this  point  as  the  origin  of  new 
coordinates;  and  (3)  to  use  the  following  form  of  the  shape  of 
the  ID  density  in  the  new  coordinates: 

N(z)  =  b  exp(-  arz2)[  1  -  ( z/z± )2], 


FIG.  5.  The  ID  spatial  photon  density  at  time  t~2llr/v,  obtained 
by  the  reshaped  form  Eq.  (33)  (solid  curve)  and  the  Gaussian  form 
(dashed  curve),  compared  with  that  of  the  Monte  Carlo  simulation 
(dots).  The  Heyney-Greenstein  phase  function  with  g=0.9  is  used, 
and  1  /la=  0.  In  the  figure,  the  unit  on  the  z  axis  is  lXr;  Rc  is  the  center 
position  of  the  distribution  computed  by  the  cumulant  solution;  zc  is 
the  distance  between  the  origin  of  the  new  coordinates  and  the 
source. 

_  j  z  >  0, 

vt  *+-  zc  i  (34) 

[z  <  0. 

At  the  ballistic  limit  z=z±,  N(z)  reduces  to  zero,  and  N(z) 
=0  when  z  is  outside  of  z±.  The  parameter  b  in  Eq.  (33)  can 
be  determined  by  normalization;  the  parameters  (a,zc)  can 
be  determined  by  fitting  the  center  and  half-width  of  the 
distribution.  This  fit  requires 

J  N(z)dz  =  1 ,  (35) 

(z}  =  J  zN(z)dz  =  Rcz-zc,  (36) 

(F2)C=J  (z—  (z))2N(z)dz  =  2Dzzvt.  (37) 

The  integrals  in  Eqs.  (35)— (37)  can  be  analytically  per¬ 
formed;  they  are  related  to  the  standard  error  function  (the 
incomplete  gamma  function,  or  the  confluent  hypergeometric 
function  of  the  first  kind): 

rP  1/2 

F{0\P)  =  e~y  dy  =  — — erf(/3),  (38) 

Jo  J- 

F{[)(0}=  f  e  y2y  dy=~[l~  , 

Jo  ^ 


where 


(33) 


(39) 
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Filn)(0)  =  e~yly1"dy=X-[(2n  -  l)F,2"~2)(/3)  -  /?2"'e^], 
Jo  2 

(40) 

F<2"+l)(/3)  =  ^  e^2y2"+'c/y  =  -[2  «F<2"_1)(/3)  -  )S2"e^]. 

J0  2 

(41) 

For  nonlinear  fitting  a  difficulty  is  how  to  quickly  find  a 
global  minimum.  The  optimization  codes  require  setting  a 
good  initial  value  of  the  parameters,  so  the  obtained  local 
minimum  is  the  true  global  minimum.  Since  we  have  no 
experience  for  setting  good  initial  parameters  at  a  special 
time,  the  following  procedure  is  used.  At  the  long  time  limit 
zc~Rcz  and  a2  =  (4 D,zvt)~\  the  distribution  approaches  the 
original  Gaussian  distribution.  We  make  a  nonlinear  fitting  at 
a  point  of  long  time  lm,  using  these  values  of  the  parameters 
as  initial  values.  Then,  we  make  a  fitting  at  -A/, 

where  At  is  a  small  time  interval,  with  the  initial  values  of 
parameters  from  those  obtained  at  tm,  to  produce  parameters 
at  Step  by  step,  the  parameters  in  the  whole  time  period 
can  be  computed.  Our  test  shows  that  the  fitting  program 
using  this  procedure  runs  quickly,  with  very  small  fitting 
error,  up  to  a  certain  short  time  limit. 

The  solid  curve  in  Fig.  5  shows  the  reshaped  spatial  dis¬ 
tribution,  Eq.  (33),  of  ID  density  at  time  t=2ltr/v,  using  the 
Heyney-Greenstein  phase  function  with  g=0.9,  which  satis¬ 
fies  causality  and  matches  the  Monte  Carlo  result  much  bet¬ 
ter  than  the  Gaussian  distribution. 

2.  3D  density 

In  this  case  the  ballistic  limit  is  represented  by  a  sphere 
with  center  located  at  the  source  position  and  the  radius  vt. 
We  move  the  peak  position  of  the  distribution  from  Rcz  to  zc 
along  the  s0=z  direction,  take  this  point  as  the  origin  of  new 
coordinates,  and  use  the  following  form  of  the  shape  of  the 
3D  density  as  a  function  of  position  in  the  new  coordinates, 
r: 

N{r)  =  b  exp[-  a(x)2r][  1  -  (F/F*)2],  (42) 

and  N{r)=0  when  r>r>,  where  \  is  the  polar  angle  of  r  in 
the  new  coordinates,  and  F*  is  the  distance  between  the  new 
origin  and  the  point  obtained  by  extrapolating  r  to  the  sur¬ 
face  of  the  ballistic  sphere, 

F*  =  [M2  -  z2  sin2(x)]1/2  -  cos(*)zc.  (43) 
In  Eq.  (42)  a(x)  is  defined  by 

«(B2  =  <4  cos2(y)  +  «±  sin2(£) .  (44) 

The  parameters  b  can  be  determined  by  normalization;  the 
parameters  ( az,a1,zc )  are  determined  by  fitting  the  center 
and  half-width  of  the  distribution.  This  fit  requires 

(z)  =  f  7  cos(x)N(r)dir  =  Rc:  -  zc,  (45) 


FIG.  6.  Time-resolved  profile  of  3D  photon  density,  where  the 
detector  is  located  at  z=3/lr  from  the  source  along  the  incident 
direction,  obtained  by  the  reshaped  form  Eq.  (42)  (solid  curve)  and 
the  Gaussian  form  (dashed  curve),  compared  with  that  of  the  Monte 
Carlo  simulation  (dots).  The  Heyney-Greenstein  phase  function 
with  g=  0.9  is  used,  and  the  absorption  coefficient  l//a=0. 

J  [F  cos(^)  -  (z)]2N(r)d3r  =  2Dzzvt,  (46) 

j  [7Mx)fmd>r-4D^.  (47) 

In  the  above  integral  dir=2'JTrzdr  d  cos(a),  integration  over  F 
can  be  analytically  performed,  and  integration  over  x  is  per¬ 
formed  numerically. 

Figure  6  shows  the  computed  time  profile  of  the  3D  den¬ 
sity  N(r,t),  with  the  source  at  the  origin  and  the  detector 
located  at  r=(0,0,3/tr),  using  the  Heyney-Greenstein  phase 
function  with  g=0.9.  The  solid  curve  is  for  the  reshaped 
form  using  Eq.  (42).  The  dashed  curve  is  for  the  Gaussian 
form,  and  the  dots  are  for  the  Monte  Carlo  simulation.  The 
results  clearly  demonstrate  an  improvement  by  use  of  the 
reshaped  form  over  use  of  the  Gaussian  form.  The  nonzero 
intensity  before  tb=3llT/v  in  the  reshaped  form  has  been 
completely  removed,  while  the  Gaussian  distribution  has 
nonzero  components  before  tb.  The  reshaped  time  profile 
matches  with  the  result  of  the  Monte  Carlo  simulation  in 
most  of  the  time  period,  but  the  peak  values  are  about  20% 
lower.  The  errors  are  much  smaller  than  those  of  the  Gauss¬ 
ian  distribution.  By  integration  over  time,  the  density  for  the 
steady  state  can  be  obtained.  The  difference  in  the  steady 
state  density  between  the  reshaped  analytical  model  and  the 
Monte  Carlo  simulation  is  about  3%. 

3.  Distribution  function  fo\r,s,t) 

When  the  detector  is  located  less  than  8/tr  from  the  source 
in  a  medium  with  large  g  factor,  the  distribution  function 
is  highly  anisotropic,  and  the  intensity  received 
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FIG.  7.  Schematic  diagram  describing  the  geometry  of  the  par¬ 
ticle  spatial  distribution  for  scattering  along  a  direction  s  A  s0.  At  a 
certain  time  t,  the  center  of  the  distribution  is  located  at  rf.  The 
half-width  of  the  spread  is  characterized  by  an  ellipsoid  (the  gray 
area).  The  large  sphere  represents  the  ballistic  limit.  The  origin  of 
the  new  coordinates  is  set  by  extending  from  |A|  to  zc.  r*  is  the 
point  obtained  by  extrapolating  a  position  r  (in  the  new  coordi¬ 
nates)  to  the  surface  of  the  ballistic  sphere,  and  the  length  r*  is 
determined  by  Eq.  (43). 


strongly  depends  on  the  angle.  One  needs  to  use  the  photon 
distribution  function  /s)(r,s,/)  instead  of  the  photon  density 
N(r,t). 

In  this  case  the  center  position  rc,  as  a  function  of  (s ,  t),  is 
not  located  on  the  axis  at  incident  direction  s0.  Without  loss 
of  generality,  we  set  the  scattering  plane  (s,s0)  as  the  x-o-z 
plane.  The  center  position  now  is  located  at  rc=(r°,0,r°). 
The  orientations  and  lengths  of  the  axes  of  the  ellipsoid, 
which  characterize  the  half-width  of  the  spread  of  the  distri¬ 
bution,  can  be  computed  as  follows.  The  nonzero  compo¬ 
nents  for  the  second  cumulant  now  are  (Bxx,Bxz,B..,By}), 
expressed  in  Eq.  (11).  Byy  represents  the  length  of  one  axis  of 
the  ellipsoid,  perpendicular  to  the  scattering  plane.  By  diago¬ 
nalizing  the  matrix 


(48) 


the  lengths  and  directions  of  the  other  two  axes  of  the  ellip¬ 
soid  on  the  scattering  plane  can  be  obtained.  In  fact,  calcu¬ 
lation  shows  that  the  direction  of  rc  is  also  the  direction  of 
one  axis  of  the  ellipsoid,  since  at  a  certain  time  t  the  direction 
rc  can  replace  s  as  the  unique  special  direction  in  the  scat¬ 
tering  plane.  In  order  to  reshape  the  distribution  we  choose  a 
new  z  axis  along  the  rc  direction,  and  move  the  peak  position 
of  the  distribution  from  1^1  to  zc,  taking  this  point  as  the 
origin  of  new  coordinates  (x,y=y,z),  as  shown  schemati¬ 
cally  in  Fig.  7. 

In  the  new  coordinates  we  use  a  shaped  form  similar  to 
that  of  the  3D  density  Eqs.  (42),  while  a{y)  in  Eq.  (42)  is 


FIG.  8.  Time-resolved  profile  of  photon  distribution  function, 
for  light  directions  0=  (a)  0  and  (b)  30°,  where  the  detector  is 
located  at  z=3/tr  from  the  source  along  the  incident  direction,  ob¬ 
tained  by  the  reshaped  form  Eq.  (42)  (solid  curves)  and  the  Gauss¬ 
ian  form  (dashed  curve),  compared  with  that  of  the  Monte  Carlo 
simulation  (dots).  The  Heyney-Greenstein  phase  function  with  g 
-0.9  is  used,  and  the  absorption  coefficient  l//„=0. 


a(x,$)2  =  <4  cos2(x>  +  a~  sin2(*)cos2($ 

+  a|sin2O0  sin2(£),  (49) 

where  x  and  <f>  are,  separately,  the  polar  angle  and  the  azi¬ 
muthal  angle  of  position  r  in  the  new  coordinates.  The  pa¬ 
rameters  (ax,ay,az,zc)  are  determined  by  fitting  the  center 
|rc|  and  lengths  of  the  three  axes  of  the  ellipsoid  character¬ 
izing  the  half-width  of  the  distribution.  In  many  cases,  the 
ellipsoid  can  be  approximately  treated  as  an  ellipsoid  of 
revolution;  the  length  of  the  axis  of  the  ellipsoid  along  the  x 
direction  is  approximately  equal  to  that  along  the  y  direction, 
and  thus  the  computation  can  be  simplified.  The  reshaped 
distribution  function  ft\r,s,t)  for  a  certain  direction  s  is 
normalized  to  F^(s,So,/). 

Figure  8  shows  the  computed  time  profile  of  the  distribu¬ 
tion  function  /j)( r,s,f),  when  the  detector  is  located  at  3/,r  in 
front  of  the  source,  using  the  Heyney-Greenstein  phase  func¬ 
tion  with  g=0.9.  Figures  8(a)  and  8(b)  are,  separately,  for 
different  directions  of  light  s:#=0  and  30°.  The  solid  curves 
are  for  the  reshaped  form  using  Eq.  (42)  and  the  dashed 
curve  is  for  the  Gaussian  form.  The  dots  are  for  the  Monte 
Carlo  simulation.  Anisotropic  distribution  is  shown  by  com¬ 
paring  Figs.  8(a)  and  8(b).  The  reshaped  distribution  removes 
intensity  before  tb=3lu/v,  which  appears  in  the  Gaussian  dis¬ 
tribution.  The  reshaped  distribution  matches  the  Monte  Carlo 
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t  (unit  of  l  Jo) 


FIG.  9.  Time-resolved  profile  of  photon  distribution  function, 
for  light  directions  6=  (a)  0  and  (b)  30°,  where  the  detector  is 
located  at  z=4/lr  from  the  source  along  the  incident  direction.  Other 
parameters  are  the  same  as  in  Fig.  8. 

result  much  better  than  the  Gaussian  distribution,  but  the 
peak  value  is  about  40%  lower  than  that  of  the  Monte  Carlo 
simulation.  Integrating  over  time  shows  that  the  difference  in 
the  steady  state  distribution  function  between  the  reshaped 
analytical  model  and  the  Monte  Carlo  simulation  is  about 
7%.  The  ratio  of  the  peak  value  at  0=30°  is  about  60%  of 
that  at  0=0,  which  shows  stronger  anisotropy  at  d=3ltr  com¬ 
pared  to  that  at  d=5lu  shown  in  Fig.  2. 

Figure  9  shows  the  distribution  function  fs\r ,s,t)  when 
the  detector  is  located  at  4/tr  in  front  of  the  source.  The 
reshaped  distribution  matches  the  Monte  Carlo  result  much 
better  than  that  at  3/tr.  It  shows  that  the  peak  intensity  at  4/tr 
is  about  one  order  of  magnitude  smaller  than  that  at  3/tr,  but 
intensity  decays  with  time  more  slowly  at  4/tr  than  at  3/tr. 

V.  DISCUSSION 

While  causality,  together  with  the  correct  center  and  half¬ 
width  of  the  distribution,  are  major  controlling  factors  in 
determining  the  shape  and  the  range  of  the  particle  distribu¬ 
tion,  the  detailed  shapes  are,  to  some  extent,  different  in  the 
different  models.  Our  choice  of  the  reshaped  form  is  based 
on  simplicity  and  ease  of  computation,  which  obviously  is 
not  the  only  available  choice.  The  initial  results  show  that  for 
g=0.9  the  parameters  in  our  model  can  be  quickly  obtained 


using  the  above  fitting  procedure  up  to  4/tr/u  for  the  3D 
case  (up  to  t^2ltr/v  for  ID  density).  The  Monte  Carlo  simu¬ 
lation  is  more  time  consuming  in  this  time  region.  This 
model  may  work  well  for  g<0.9  in  the  above  time  region, 
because  there  is  less  forward  scattering  for  a  smaller  g  factor. 
The  fitting  error  begins  to  increase  during  3llr/v<t<4llr/v. 
At  early  time  t<3ljv,  rc  becomes  very  close  to  the  ballistic 
limit  i )t\  the  front  edge  of  the  distribution  almost  perpendicu¬ 
larly  jumps  at  the  position  vt.  In  this  case,  the  parameter  zc 
*=vt  in  our  model,  is  difficult  to  adjust  through  the  fitting 
program.  A  more  suitable  model  in  this  early  time  period  is 
needed.  Of  cause,  Monte  Carlo  simulation  also  runs  fast  for 
short  times  and  small  spatial  regions.  For  s  at  the  near  back- 
scattering  direction,  the  Gaussian  distribution  can  be  a  good 
approximation  as  shown  in  Fig.  3,  because  most  particles 
suffer  many  scattering  events  to  transfer  from  the  forward 
direction  to  the  backward  direction.  Our  calculation  shows 
that  the  center  position  rc  is  close  to  the  source  for  6 
=  1 80°  and  far  from  the  ballistic  limit;  hence,  reshaping  has 
little  effect  on  the  backscattering  case. 

In  addition  to  improving  convergence,  separating  the  bal¬ 
listic  component  from  the  scattered  component  also  provides 
a  more  appropriate  time-resolved  profile  for  transmission.  In 
the  time-resolved  transmission  profile  the  ballistic  compo¬ 
nent  is  described  by  a  sharp  jump  exactly  at  time  vt,  sepa¬ 
rated  from  the  later  scattered  component.  The  intensity  of  the 
ballistic  component,  compared  to  the  scattered  component, 
strongly  depends  on  the  g  factor.  For  g=  0,  /tr=4,  the  ballistic 
component  decays  to  exp(-l)=0.368  at  distance  l/tr.  But  for 
g=0.9  it  decays  to  exp(-10)=4.54X  10~5  at  l/tr,  because  /tr 
=  104.  The  jump  of  the  ballistic  component  can  be  seen  in 
experiments  of  transmission  of  light  for  a  medium  of  small 
sized  scatterers  (small  g  factor),  but  is  difficult  to  observe  for 
a  medium  of  large  sized  scatterers  (large  g  factor).  Our  for¬ 
mula  presented  in  Sec.  Ill  provides  a  good  estimation  for 
both  small  and  large  g  factors  by  explicitly  separating  these 
two  components. 

Using  the  obtained  analytical  expressions,  the  distribution 
I(r,s,t)  can  be  computed  very  fast.  The  cumulant  solution  of 
the  BTE  has  been  extended  to  the  case  of  a  polarized  photon 
distribution  [11],  and  to  semi-infinite  and  slab  geometries 

[12] .  Using  a  perturbation  method,  the  distribution  I(r,s,t) 
in  a  weak  heterogeneous  medium  can  be  calculated  based  on 
the  cumulant  solution  of  the  BTE  [12].  The  nonlinear  effect 
for  strongly  heterogeneous  objects  inside  a  medium  can  also 
be  calculated  using  a  correction  of  the  “self-energy”  diagram 

[13] .  Hence,  the  analytical  form  of  the  solution  of  the  BTE 
may  have  many  different  applications. 

In  summary,  the  analytical  cumulant  solution  of  the  Bolt¬ 
zmann  transport  equation  is  improved  in  two  respects.  The 
ballistic  component  is  separated  and  the  cumulants  for  the 
scattered  component  are  computed.  This  treatment  ensures 
that  summation-over  /  is  convergent.  We  replace  the  Gauss¬ 
ian  distribution  by  a  different  shaped  form,  which  satisfies 
causality,  and  maintains  the  correct  center  position  and  the 
correct  half-width  of  the  distribution  computed  by  our  ana¬ 
lytical  formula.  Our  results  show  that  the  reshaped  distribu¬ 
tion  matches  that  obtained  by  the  Monte  Carlo  simulation 
much  better  than  does  the  Gaussian  distribution. 
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APPENDIX  2 


Electric  field  Monte  Carlo  simulation  of 
polarized  light  propagation  in  turbid 

media 

Min  Xu 

Institute  for  Ultrafast  Spectroscopy  and  Lasers,  New  York  State  Center  of  Advanced 
Technology  for  Ultrafast  Photonics,  and  Department  of  Physics,  The  City  College  and 
Graduate  Center  of  City  University  of  New  York,  New  York,  NY  10031 


Abstract:  A  Monte  Carlo  method  based  on  tracing  the  multiply  scattered 
electric  field  is  presented  to  simulate  the  propagation  of  polarized  light  in 
turbid  media.  Multiple  scattering  of  light  comprises  a  series  of  updates  of 
the  parallel  and  perpendicular  components  of  the  complex  electric  field 
with  respect  to  the  scattering  plane  by  the  amplitude  scattering  matrix  and 
rotations  of  the  local  coordinate  system  spanned  by  the  unit  vectors  in  the 
directions  of  the  parallel  and  perpendicular  electric  field  components  and 
the  propagation  direction  of  light.  The  backscattering  speckle  pattern  and 
the  backscattering  Mueller  matrix  of  an  aqueous  suspension  of  polystyrene 
spheres  in  a  slab  geometry  are  computed  using  this  Electric  Field  Monte 
Carlo  (EMC)  method.  An  efficient  algorithm  computing  the  Mueller  matrix 
in  the  pure  backscattering  direction  is  detailed  in  the  paper. 
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1.  Introduction 

The  propagation  of  polarized  light  in  turbid  media  is  fundamental  to  many  practical  appli¬ 
cations  of  considerable  interest  including  remote  sensing  of  clouds  and  imaging  of  colloidal 
suspensions  and  biological  materials]!^!].  Due  to  the  lack  of  analytical  solutions  to  radiative 
transfer[5]  of  polarized  light  within  a  bounded  medium,  numerical  solutions[6-8]  of  the  radia¬ 
tive  transfer  equation  and  Monte  Carlo  simulations[9-18]  of  propagation  of  polarized  light  in 
turbid  media  have  been  pursued  extensively  and  applied  to  characterization  of  clouds,  particle 
suspensions,  and  biological  materials. 

Most  Monte  Carlo  methods  trace  the  Stokes  vector  I  =  ( I,Q,U,V)T  in  simulation  where 
I  =  (|£,|2  +  \Er\ 2),  Q  =  (|£/|2  -  |£r|2},  U  =  (EJEr+E,E*),  V  =  -i (EJ Er  - E,E;),  E,  and 
Er  are  two  orthogonal  complex  electric  field  components  perpendicular  to  the  propagation  di¬ 
rection,  the  superscript  “7”’  denotes  transpose,  and  {)  means  ensemble  average.  Light  scat¬ 
tering  involves  a  rotation  of  the  Stokes  vector  to  a  local  scattering  reference  frame  and  the 
multiplication  of  the  Stokes  vector  by  the  4x4  Mueller  matrix  which  prescribes  how  polar¬ 
ized  light  is  scattered  by  an  isolated  particle  in  that  reference  frame.  The  Stokes  vector  is 
summed  up  at  the  detector  assuming  the  detected  light  is  incoherent.  Many  implementations 
of  the  above  Monte  Carlo  approach  to  simulate  polarized  light  propagation  in  turbid  media 
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have  appeared  in  the  literature[9-18].  Some  recent  work  includes  [13, 14, 16-18],  Symmetry 
considerations^  9],  randomly  polarized  incident  light[  1 7],  next  event  point  estimators[17]  and 
other  techniques[14, 16]  have  been  investigated  to  improve  the  efficiency  of  the  Monte  Carlo 
computation.  Various  implementations  of  the  above  Monte  Carlo  approach  have  also  been  com¬ 
pared^]. 

In  this  article,  we  present  a  Monte  Carlo  method  based  on  tracing  the  multiply  scattered 
electric  field  to  simulate  the  propagation  of  polarized  light  in  turbid  media.  Multiple  scattering 
of  light  comprises  a  series  of  updates  of  the  parallel  and  perpendicular  components  of  the 
complex  electric  field  with  respect  to  the  scattering  plane  by  the  amplitude  scattering  matrix 
and  rotations  of  the  local  coordinate  system  spanned  by  the  unit  vectors  in  the  directions  of 
the  parallel  and  perpendicular  electric  field  components  and  the  propagation  direction  of  light. 
The  phase  of  light  accrues  from  phase  delays  due  to  the  interaction  of  light  with  scatterers  and 
propagation  through  the  host  medium.  In  contrast  with  the  conventional  Monte  Carlo  approach 
based  on  tracing  the  Stokes  vector,  the  Electric  Field  Monte  Carlo  (EMC)  method  traces  the 
electric  field  and  phase  of  light  and  makes  it  possible  to  simulate  also  coherent  properties  of 
multiply  scattered  light.  The  algorithm  of  EMC  is  straightforward  and  can  be  easily  adapted  to 
simulate  the  propagation  of  polarized  light  in  optically  active  medium. 

After  outlining  the  theoretical  formalism  of  the  Electric  Field  Monte  Carlo  method  in  sec¬ 
tion  2,  we  present  EMC  computations  of  the  backscattering  speckle  pattern  and  the  backscat- 
tering  Mueller  matrix  of  an  aqueous  suspension  of  polystyrene  spheres  in  a  slab  geometry  in 
section  3.  Special  considerations  to  improve  the  efficiency  of  computing  the  Mueller  matrix  in 
the  pure  backscattering  direction  are  detailed  in  section  3.  We  conclude  in  section  4. 

2.  Theoretical  formalism 

Light  scattering  by  small  particles  is  succinctly  described  by  the  amplitude  scattering  ma- 
trix[2 1  ].  For  spherical  or  randomly-oriented  aspheric  scatterers,  the  amplitude  scattering  matrix 
is  diagonal  and  depends  only  on  the  scattering  angle  6  due  to  the  symmetry [22].  The  parallel 
and  perpendicular  components  of  the  electric  field  with  respect  to  the  scattering  plane  are  scat¬ 
tered  according  to  Sj(0)  where  j  —  2,1,  respectively[21]. 

The  scattering  of  photons  takes  a  simple  form  in  the  local  orthonormal  coordinate  system 
(m.n,s)  where  m  and  n  are  the  unit  vectors  in  the  directions  of  the  parallel  and  perpendicular 
components  of  the  electric  field  with  respect  to  the  scattering  plane  of  the  previous  scattering 
event  and  s  is  the  photon’s  propagation  direction  prior  to  the  current  scattering  [see  Fig.  1],  The 
propagation  direction  s'  of  the  photon  after  the  current  scattering  is  given  by: 

s'  =  msin0cos<|>  4-  nsin  Osin  <p  +  scos  0  (1) 

where  6  and  </>  are  the  scattering  and  azimuthal  angles  of  the  current  scattering  respectively. 

The  current  scattering  plane  is  the  plane  spanned  by  s  and  s' .  The  unit  vectors  in  the  direc¬ 
tions  of  the  parallel  and  perpendicular  electric  fields  with  respect  to  the  current  scattering  plane 
are  given  by 


ei  =  mcos0  +  nsm<p  (2) 

e2  =  —  m  sin  0  +  n  cos  0 


prior  to  scattering  and 


o',  =  ji/mcos<]>  +  jifnsin0 -sinOs, 


(3) 
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Fig.  1 .  A  photon  moving  along  s  is  scattered  to  s'  with  a  scattering  angle  0  and  an  azimuthal 
angle  0  inside  a  local  coordinate  system  spanned  by  orthonormal  bases  (m,n,s).  ei .2  and 
e'j  1  are  the  unit  vectors  parallel  and  perpendicular  to  the  current  scattering  plane  spanned 
by  s  and  s'  prior  to  and  after  scattering.  The  local  coordinate  system  (m.n.s)  is  rotated  to 
(in',  n',  s')  after  scattering. 


after  scattering,  respectively.  The  local  coordinate  system  (m.n.s)  is  rotated  to  (m',n',s') 
whose  m'  =  e',  and  n'  =  e'2  after  scattering.  The  incident  electric  field  E  =  £|tn  +  £2 n  is 
scattered  to  E'  =  E[  rn'  +  E'2n'  whose  parallel  and  perpendicular  components  are  given  by 
E\  —  S2E  •  ei  and  E'2  =  St  E  •  e2,  respectively. 

We  can  now  summarize  the  updating  rules  of  the  local  coordinate  system  and  the  electric 
field  in  EMC.  For  each  scattering  with  the  scattering  angle  0  and  the  azimuthal  angle  <p,  the 
local  coordinate  system  is  updated  by 


with 


A  = 


cos  0  COS  0 
— sin$ 
sin  0  cos  </> 


cos  0  sin  § 
cos  <p 
sin  0  sin  0 


and  the  electric  field  by 


(5) 

(6) 
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with 


(7) 


8=[F(e.*r,/2(  %cosi  ?sint  V 

L  v  n  —Sj  sin  0  Si  cos  0  J 

We  have  introduced  an  extra  factor  F(0, 0),  the  scattered  wave  intensity  propagating  along  the 
(0,0)  direction,  given  by 

£(0,0)  =  (|S2|2 cos2  0  +  |Si |2 sin2 0)  |£,  |2  +  (|S2 12 sin2 0  +  |S,  |2 cos2  0)  |£2|2 

+2  (|S2|2  -  |S,|2)  cos 0  sin 091  [£,  (£2)*]  (8) 

to  normalize  the  light  intensity.  The  intensity  of  the  incident  light  \E\  |2  + |£2|2  =  1  is  assumed 
unity  in  Eq.  (8).  The  scattered  light  intensity  j£[|2  +  \E'2\2  is  then  conserved  (equal  to  unity) 
in  the  series  of  scattering  events.  Absorption  of  light,  if  any,  is  accounted  for  by  adjusting  the 
photon  weight  as  in  the  conventional  Monte  Carlo  simulations[23]. 

The  orthogonal  matrix  A  rotates  the  local  coordinate  system  (m,n,s)  each  time  the  pho¬ 
ton  is  scattered  by  a  scatterer.  The  electric  field  is  updated  simultaneously  by  the  matrix  B. 
The  free  path  between  consecutive  scattering  events  is  sampled  in  the  same  fashion  as  that  in 
a  conventional  Monte  Carlo  method.  The  state  of  a  multiply  scattered  photon  is  specified  by 
the  final  local  coordinate  system  ^mW,nW,sWj  after  consecutive  rotations,  the  final  com¬ 
plex  electric  field  components  E\n\  after  consecutive  updates,  and  the  optical  path  /  inside  the 
host  medium  where  n  denotes  the  number  of  scattering  events  the  photon  have  experienced  be¬ 
fore  being  detected.  The  detected  electric  field  is  given  by  E</  =  exp (ikl) 

where  k  =  2n/X  is  the  wave  number,  A  is  the  wavelength  of  light  in  the  host  medium,  and  we 
have  assumed  a  convention  that  the  temporal  dependence  of  a  plane  wave  of  frequency  (0  is 
exp  (-/to/).  The  phase  of  the  detected  photon  accrues  from  phase  delays  due  to  the  interaction 
of  light  with  scatterers  and  propagation  through  the  host  medium.  The  photon  weight  w,  unity  at 
incidence,  is  multiplied  by  the  albedo  of  the  scatterer  at  each  scattering  event.  Once  the  photon 
hits  a  detector,  the  electric  field  at  the  detector  is  increased  by  w’^E c/  and  the  Stokes  vector  is 
increased  by  the  Stokes  vector  computed  from  the  electric  field  w1/2 E(/. 

Propagation  of  multiply  scattered  light  is  describable  by  the  radiative  transport  equation  in 
which  any  interference  effects  are  neglected[5].  The  conventional  Monte  Carlo  methods  based 
on  tracing  the  Stokes  vector  of  light  solve  the  radiative  transfer  equation  numerically  and  do  not 
include  any  correlation  effect  of  multiply  scattered  light.  By  tracing  the  electric  field  associated 
with  a  wave  packet,  EMC  provides  much  more  detailed  information  about  the  propagation  of 
multiply  scattered  light  than  just  the  transfer  of  energy.  EMC  should  be  regarded  as  a  method 
to  sample  the  probability  distribution  of  the  electric  field  at  selected  spatial  positions  where 
the  detectors  are  located.  Detectors  of  finer  resolution  than  the  speckle  size  are  required  to 
resolve  the  interference  pattern  of  light  well.  Detectors  of  a  larger  cell  size  will  smear  the 
interference  pattern.  It  should  be  noted  that  EMC  only  probes  one  instantaneous  picture  of  the 
disordered  medium  when  accumulating  the  electric  field  coherently  while  it  yields  the  detected 
light  ensemble  averaged  over  all  pictures  of  the  disordered  medium  when  accumulating  the 
Stokes  vector  incoherently.  This  point  is  addressed  in  more  detail  in  section  3.1. 

One  key  step  in  the  Monte  Carlo  simulation  of  polarized  light  is  the  sampling  of  the  scat¬ 
tering  angles  (0,0)  which  distribute  according  to  the  normalized  phase  function  p(6,<p)  = 
F(6.<t>)/KQsalx 2  where  gSC3  is  the  scattering  efficiency,  x  =  ka  is  the  size  parameter,  and  a 
is  the  radius  of  the  particle.  The  rejection  technique[24, 25]  has  been  widely  used  to  sample 
such  a  distribution  function.  The  procedure  is  to  choose  a  doublet  (ju  =  cos  0,0)  where  p  and 
0  are  uniformly  distributed  over  [—1,1]  and  [0,2/t:]  respectively  and  a  random  number  /  uni- 
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fonnly  distributed  over  [O.maxg  ^ p(6.  0)],  the  doublet  (p,  0)  is  accepted  as  the  scattering  and 
azimuthal  angles  if /  <  /?(arccosjti,0),  otherwise  generate  a  new  doublet  {p,<p),  a  new  /  and 
start  over.  Each  trial  involves  one  Mie  calculation  of  the  amplitude  scattering  matrix  at  that  trial 
scattering  angle.  The  number  of  rejections  per  scattering  event  is  large  (on  the  order  of  one 
for  Rayleigh  particles  and  can  be  much  more  significant  for  larger  particles).  This  limits  the 
efficiency  of  the  Monte  Carlo  simulation.  Mie  calculations  of  the  amplitude  scattering  matrix 
are  hence  usually  performed  on  a  fixed  set  of  scattering  angles  in  advance  to  generate  a  table  of 
scattering  matrices  and  reduce  the  computation  load. 

We  generate  the  scattering  angle  0  by  drawing  of  one  random  number  uniformly  distributed 
over  [0. 1]  and  looking  up  an  inverse  table  of  the  marginal  distribution  function 


P{0)=  [  p{8,<l>)d<l>  = 
Jo 


l&P  +  M2 

()scn-V“ 


(9) 


which  does  not  depend  on  the  electric  field  and  is  pre-calculated  before  simulation.  The  az¬ 
imuthal  angle  0  is  then  obtained  by  the  rejection  method  for  the  conditional  probability 
/;(<;') 1 0)  =  p(0 ,  <p) / p(6) .  One  Mie  calculation  for  the  looked  up  scattering  angle  is  required 
if  the  table  of  scattering  matrices  is  not  generated  in  advance  for  each  scattering  event. 

A  different  sampling  strategy  is  to  sample  the  scattering  angle  0  according  to  p(0)  while 
to  sample  0  uniformly  distributed  over  [0, 2tt].  At  the  same  time,  we  modify  F(0,0)  in  the 

update  rule  of  the  electric  field  (6)  to  F'(6)  =  +  |Si  | 2  j  /2  such  that  the  light  intensity  is 

no  longer  conserved  in  the  series  of  scattering  events.  This  strategy  saves  the  rejection  sampling 
of  the  angle  0  but  the  simulation  result  is  prone  to  be  contaminated  by  hotspots  and  has  a  larger 
variance  compared  to  the  first  strategy  because  it  is  unfavorable  to  the  statistical  variance  when 
photons  of  varying  weights,  rather  than  equal  weights,  are  accumulated  by  the  detector.  The 
second  sampling  strategy  is  not  used  in  simulation. 


3.  Results  and  Discussion 

3.1.  Backscattering  speckle  pattern 

Due  to  the  large  ratio  between  the  velocities  of  light  and  of  the  scatterers,  coherent  light 
probes  an  instantaneous  picture  of  the  disordered  medium  and  realizes  speckle  patterns  when 
the  multiple  scattered  light  emerges  from  within  the  medium.  The  a  =x ,y,z  component  of 
the  outgoing  light  intensity  7a(0,0)  in  direction  (0,0)  on  the  boundary  of  the  medium  com¬ 
prises  a  multitude  of  independently-phased  additive  complex  electric  fields  in  that  direction 
(outside  the  coherent  backscattering  cone[26-28])  and  follows  a  Rayleigh  distribution[29] 
p(Ia)  —  exp(-/„/  (/a})  where  (I,,)  is  the  mean  intensity.  The  nonualized  light  intensity 

7]  =  laj  {/„)  follows  the  exponential  distribution  exp(— r\). 

We  perform  EMC  to  study  the  plane  wave  light  backscattering  from  an  aqueous  solution  of 
polystyrene  spheres  inside  a  slab.  Incident  light  is  polarized  in  the  x  direction  and  propagates  in 
the  z  direction  (normal  to  the  surface  of  the  slab).  The  size  of  the  polystyrene  sphere  is  0.46/tm. 
The  wavelength  of  the  incident  wave  is  2,vac  =  0.52p  m  in  vaccum.  The  refractive  indices  of  the 
polystyrene  sphere  and  water  is  «sct  =  1.59  and  =  1 .33,  respectively.  The  mean  scattering 
length  of  light  inside  the  solution  is  ls  =  2.80/im  with  2nn^ls/Xy.ic  =  45.  The  thickness  of 
the  slab  is  Lz  =  4 !s.  Total  5  x  108  photons  are  launched  into  the  medium.  Both  backscattering 
electric  field  and  Stokes  vector  are  recorded  versus  the  direction  (0,0)  of  the  backscattered 
light  where  0  is  the  angle  away  from  the  surface  normal  (0  <  0  <  n/2)  and  0  is  the  azimuthal 
angle  (0  <  0  <  2k).  The  ranges  of  the  angles  0  and  0  are  uniformly  divided  into  1 125  and  360 
bins  in  simulation,  respectively. 
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(a) 


(b) 


Fig.  2.  (a)  Speckle  pattern  formed  by  the  angular-resolved  backscattering  light,  (b)  Nor¬ 
malized  speckle  4/ (4)  follows  a  negative  exponential  distribution. 


EMC,  like  an  experiment  probing  one  static  random  medium,  only  probes  one  instanta¬ 
neous  picture  of  the  disordered  medium.  The  ensemble  average  is  realized  in  experiments,  for 
example,  through  the  movement  of  the  scatterers  in  the  medium  or  sampling  different  regions 
of  the  random  media.  EMC,  as  a  numerical  experiment,  can  record  both  the  electric  field  and 
the  Stokes  vectors  in  simulation.  The  recorded  electric  field  added  coherently  yields  the  emer¬ 
gent  light  from  the  one  instantaneous  picture  of  the  disordered  medium.  The  recorded  Stokes 
vector  added  incoherently,  on  the  other  hand,  yields  the  emergent  light  ensemble  averaged  over 
all  pictures  of  the  disordered  medium. 

Figure  2(a)  displays  the  x  component  /*/(/*)  where  4  =  \EX\2  and  (Ix)  is  estimated  by 
the  mean  of  the  first  and  second  Stokes  vectors.  The  appearance  of  speckles  is  obvious  in 
Fig.  2(a).  The  first  order  statistics  about  4/  (4)  can  be  computed  from  its  histogram.  This 
yields  a  negative  exponential  distribution  exp (— rj)  as  expected  [see  Fig.  2(b)].  The  other  twoy 
and  z  components  of  light  behave  similarly  and  not  displayed. 

3.2.  Backscattering  Mueller  matrix 

We  then  compute  the  backscattering  Mueller  matrix  from  a  pencil  beam  incident  normally  on 
an  aqueous  solution  of  polystyrene  spheres  inside  a  slab.  The  size  of  the  polystyrene  sphere  is 
2jum.  The  wavelength  of  the  incident  wave  is  Avac  =  0.63/ini  in  vacuum.  The  refractive  indices 
of  the  polystyrene  sphere  and  water  is  «sci  =  1.59  and  «t,g  =  1.33,  respectively.  The  thickness 
of  the  slab  is  Lz  =  44-  Total  3  x  108  photons  are  launched  into  the  medium. 

To  improve  its  efficiency  of  the  Monte  Carlo  simulation,  we  combined  symmetry  consid¬ 
erations^],  randomly  polarized  incident  light[17]  and  the  next  event  point  estimator[  1 7]  in 
computing  the  backscattering  Mueller  matrix.  The  computation  time  is  less  than  2  hours  for 
each  108  photons  launched  using  one  Pentium  111  1GHz  CPU.  In  computation  of  the  Stokes 
vector  in  Monte  Carlo  simulations,  one  should  note  that  the  Stokes  vector  upon  which  Mueller 
matrix  is  defined  depends  on  the  reference  system  used.  The  backscattered  Stokes  vector  is 
defined  in  the  reference  system  whose  x'yJz  axes  coincide  with  —x.y,  — z  axes  of  the  reference 


#5731  -S15.00US 
(C)  2004  OSA 


Received  1 5  November  2004;  revised  1 3  December  2004;  accepted  1 4  December  2004 
27  December  2004 /Vol.  12,  No.  26  /  OPTICS  EXPRESS  6536 


system  of  the  incident  Stokes  vector  respectively. 

In  the  backscattering  geometry  investigated  here  where  the  incident  light  is  normal  to  the 
surface  of  the  medium  and  the  outgoing  beam  is  exactly  backscattering,  the  Mueller  matrix 
Afbs(p,  p)  of  the  medium  shall  satisfy  the  following  relationf  1 9, 30] 

iri°s{p4 >)  =  /?(— p)A/bs(p,  p  =  o)/?(-p)  (io) 


where  p,  p  is  the  polar  coordinate  of  the  position  of  the  outgoing  beam  and  R(<p)  is  the  rotation 
matrix 
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The  reduced  Mueller  matrix  Mqs(p)  =  AT^p.  p  =  0)  relates 


1'  =  A^s(p)l' 


(12) 


where  I'  0  are  the  incident  and  outgoing  Stokes  vectors  with  respect  to  the  p  =  0  plane. 

In  our  simulation,  the  incident  electric  field  is  randomly  polarized  as  E,  =  cos  ae~'^x  + 
sin  ae'Py  where  a  €  (0.  n/2)  and  |3  £  (0,  K)  are  unifonn  random  numbers.  The  corre¬ 
sponding  incident  Stokes  vector  with  respect  to  the  xy  axes  (i.e.,  the  y  =  0  plane)  is  I,-  = 
(l,cos2a,sin2acos2j3,sin2asin2/3)T .  The  incident  Stokes  vector  with  respect  to  the  p  =  0 
plane  is  given  by  I-  =  /?(p)I,-.  The  ensemble  average  of  Ij  (I')7  over  the  random  variates  a  and 
1 3  yields 


[;(i;)r)  =  /j(P)(i,if)7?(-P)  =  /?(P) 
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Noting  that  the  inverse  of  the  matrix 


(I'iflf)  is  given  by 
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we  obtain 

<(pH1«(ID7£  (15) 

from  Eq.  (12).  The  backscattering  Mueller  matrix  is  then  computed  from  the  reduced  Mueller 
matrix  by  Eq.  (10). 

The  computed  backscattering  Mueller  matrix  is  displayed  in  Fig.  (3).  Each  matrix  element 
is  given  as  a  two-dimensional  image  of  the  surface,  20/s  x  20/,  in  size,  with  the  laser  being 
incident  in  the  center.  The  displayed  Mueller  matrix  has  been  normalized  by  the  maximum  light 
intensity  of  the  (1,1)  element  of  Mbs.  Only  7  elements  of  the  Mueller  matrix  are  independent. 
The  symmetrical  relation  between  different  elements  of  the  Mueller  matrix  has  been  considered 
by  Kattawa  et.  al.[19, 31],  The  symmetry  of  the  backscattering  Mueller  matrix  in  Fig.  (3)  agrees 
with  Kattawa  et.  al.[13, 19,31]  and  Berezhnyy  et.  al.[32].  The  elements  and  vanish  as 
predicted  by  the  theory[13, 19], 
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Fig.  3.  Backscattering  Mueller  matrix  for  the  slab.  All  4  x  4  matrix  element  are  displayed 
as  a  two-dimensional  image  of  the  surface,  20/,  x  20/,  in  size,  with  the  laser  being  incident 
in  the  center.  The  displayed  Mueller  matrix  has  been  normalized  by  the  maximum  light 
intensity  of  the  (1,1)  element.  The  parameters  of  the  slab  is  given  in  the  text. 


It  is  more  instructive,  though,  to  display  the  reduced  Mueller  matrix  M§s(p).  Each  element 
of  the  reduced  Mueller  matrix  is  given  as  a  one-dimensional  curve  versus  the  distance  p  from 
the  origin  in  Fig.  (4).  The  reduced  backscattering  Mueller  matrix  is  found  to  be  2  x  2  block 
diagonal.  The  first  two  elements  of  the  Stokes  vector  1\  2  are  then  decoupled  from  the  rest  two 
elements  of  the  Stokes  vector  I'i4.  This  means,  for  example,  the  backscattered  light  due  to  a 
normally  incident  light  linearly  polarized  in  the  0  =  0  plane  (Ij  =  (1, 1,0, 0)r)  has  no  circular 
polarization  component  (the  fourth  element  of  Stokes  vector  1'0  =  Mqs(p)IJ  is  always  zero)  with 
respect  to  the  0=0  plane.  The  circular  polarization  component  of  this  backscattered  light  is  no 
longer  zero  with  respect  to  a  different  reference  frame  rather  than  the  0  =  0  plane. 

4.  Conclusion 

In  conclusion,  we  have  presented  a  Monte  Carlo  method  based  on  tracing  the  electric  field  for 
simulating  polarized  light  propagation  in  turbid  media.  The  Electric  Field  Monte  Carlo  method 
has  been  successfully  used  to  study  the  backscattering  speckle  pattern  and  the  backscattering 
Mueller  matrix  of  an  aqueous  suspension  of  polystyrene  spheres  in  a  slab  geometry.  The  tracing 
of  the  electric  field  in  simulation  makes  the  Electric  Field  Monte  Carlo  method  possible  to 
simulate  also  coherent  properties  of  light.  The  EMC  source  code  will  be  put  in  the  public 
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Fig.  4.  Reduced  backscattering  Mueller  matrix  for  the  slab.  All  4  x  4  elements  of  the  re¬ 
duced  Mueller  matrix  is  displayed  as  a  one-dimensional  curve  versus  the  distance  p/ls 
from  the  origin.  The  reduced  backscattering  Mueller  matrix  is  2  x  2  block  diagonal. 

domain[33]. 
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APPENDIX  3 


Three-dimensional  localization  and  optical 
imaging  of  objects  in  turbid  media  with 
independent  component  analysis 


M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano 


A  new  approach  for  optical  imaging  and  localization  of  objects  in  turbid  media  that  makes  use  of  the 
independent  component  analysis  (ICA)  from  information  theory  is  demonstrated.  Experimental  arrange¬ 
ment  realizes  a  multisource  illumination  of  a  turbid  medium  with  embedded  objects  and  a  multidetector 
acquisition  of  transmitted  light  on  the  medium  boundary.  The  resulting  spatial  diversity  and  multiple 
angular  observations  provide  robust  data  for  three-dimensional  localization  and  characterization  of 
absorbing  and  scattering  inhomogeneities  embedded  in  a  turbid  medium.  ICA  of  the  perturbations  in  the 
spatial  intensity  distribution  on  the  medium  boundary  sorts  out  the  embedded  objects,  and  their  locations 
are  obtained  from  Green’s  function  analysis  based  on  any  appropriate  light  propagation  model.  Imaging 
experiments  were  carried  out  on  two  highly  scattering  samples  of  thickness  approximately  50  times  the 
transport  mean-free  path  of  the  respective  medium.  One  turbid  medium  had  two  embedded  absorptive 
objects,  and  the  other  had  four  scattering  objects.  An  independent  component  separation  of  the  signal,  in 
conjunction  with  diffusive  photon  migration  theory,  was  used  to  locate  the  embedded  inhomogeneities.  In 
both  cases,  improved  lateral  and  axial  localizations  of  the  objects  over  the  result  obtained  by  use  of 
common  photon  migration  reconstruction  algorithms  were  achieved.  The  approach  is  applicable  to 
different  medium  geometries,  can  be  used  with  any  suitable  photon  propagation  model,  and  is  amenable 
to  near-real-time  imaging  applications.  ©  2005  Optical  Society  of  America 
OCIS  codes:  170.7050,  170.5280,  170.3660,  100.3190,  100.3010. 


1.  Introduction 

Optical  tomographic  imaging  of  objects  in  turbid  me¬ 
dia  is  an  aggressively  pursued  area  of  contemporary 
research  that  derives  impetus  from  a  variety  of  po¬ 
tential  practical  applications.1-17  Of  particular  inter¬ 
est  are  medical  applications  in  which  optical 
tomography  and  spectroscopy  have  the  potential  to 
provide  diagnostic  information  about  tumors  in 
breast  and  prostate  tissue  and  functional  information 
about  brain  activities.  Simultaneous  developments  in 
experimental  apparatus  and  techniques  for  object  in¬ 
terrogation  and  signal  acquisition,1 2,4,5,18,19  analytical 
models  for  light  propagation,10-20-22  and  computer  al- 


The  authors  are  with  the  Department  of  Physics,  the  Institute 
for  Ultrafast  Spectroscopy  and  Lasers,  New  York  State  Center  of 
Advanced  Technology  for  Ultrafast  Photonic  Materials  and  Appli¬ 
cations,  City  College  of  New  York,  136th  Street  Convent  Avenue 
J419,  New  York,  New  York  10031. 

Received  14  July  2004;  revised  manuscript  received  24  Novem¬ 
ber  2004;  accepted  5  December  2004. 
0003-6935/05/101889-09$15.00/0 
©  2005  Optical  Society  of  America 


gorithms  for  image  reconstruction6-8  hold  promise  for 
realization  of  these  potentials  of  optical  tomography. 

Researchers  today  use  continuous-wave,  amplitude- 
modulated,  or  ultrashort  light  pulses  to  probe  the 
target(s)  embedded  in  the  turbid  medium  and  obtain 
steady-state,  frequency-domain,  or  time-varying  op¬ 
tical  signals,  respectively,  by  using  a  variety  of  detec¬ 
tion  schemes.2-4-5-18-19  Multiple  scattering  of  light  in 
turbid  media,  such  as  breast  tissue,  precludes  direct 
imaging  of  embedded  targets.  One  then  resorts  to  an 
inverse  image  reconstruction  (HR)6-8  approach  that 
uses  a  forward  model  for  light  propagation,  the  mea¬ 
sured  light  intensity  distribution  on  the  boundary  of 
the  turbid  medium,  and  an  inversion  algorithm  to 
generate  a  map  of  the  optical  properties,  such  as  the 
absorption  coefficient  (|xa)  and  the  scattering  coeffi¬ 
cient  (|xs),  of  the  medium  and  the  embedded  objects. 
The  objects  are  desired  to  appear  as  localized  inho¬ 
mogeneities  in  the  spatial  distribution  of  these  opti¬ 
cal  coefficients. 

The  inversion  problem  is  ill  posed,  and  its  adequate 
theoretical  treatment  is  critical  to  the  achievement  of 
a  unique  solution.8  Although  three  general 
approaches — Radon-transform-type  straight  line  in- 
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tegrals,23  modeling  of  light  scattering  as  a  Markov 
random  process,24'25  and  development  and  inversions 
of  a  partial  differential  equation  (PDE)  of  diffusion 
type — are  pursued,  it  is  the  PDE-based  methods  that 
seem  more  practical  in  consideration  of  the  signal-to- 
noise  (S/N)  ratio  of  the  data  and  computationally 
efficient  methods  available  for  solution.6  8  The  com¬ 
monly  used  PDE  is  the  diffusion  approximation  (DA) 
of  the  radiative  transfer  equation  (RTE).  Both  itera¬ 
tive  reconstruction  and  noniterative  linearized  inver¬ 
sion  approaches  have  been  used  to  solve  the  inversion 
problem,  which  is  weakly  nonlinear  with  limited  suc¬ 
cess.  The  reconstruction  of  images  with  adequate 
spatial  resolution  and  optical  contrast  and  the  deter¬ 
mination  of  the  location  of  the  inhomogeneities  re¬ 
main  formidable  tasks.  The  time  required  for  data 
acquisition  and  image  reconstruction  is  another  im¬ 
portant  consideration. 

In  this  article  we  present  a  simple  and  fast  ap¬ 
proach  that  employs  (for  test  of  concept)  continuous- 
wave  transillumination  measurements  and  a  novel 
algorithm  based  on  independent  component  analysis 
(ICA)26’27  from  information  theory  to  locate  tumorlike 
inhomogeneities  embedded  in  breast-simulating  tur¬ 
bid  media.  ICA  has  been  successfully  applied  in  a 
variety  of  other  applications.27-30  We  refer  to  this 
information-theory-inspired  approach  as  optical 
tomography  that  uses  independent  component  anal¬ 
ysis,  abbreviated  as  OPTICA.  Experimental  arrange¬ 
ment  for  OPTICA  realizes  a  multisource  illumination 
and  multidetector  signal-acquisition  scheme  that 
provides  a  variety  of  spatial  and  angular  views  that 
are  essential  for  three-dimensional  (3-D)  object  local¬ 
ization.  Multisource  illumination  is  realized  in  prac¬ 
tice  by  scanning  the  input  surface  (or  source  plane) 
across  the  incident  beam  in  a  two-dimensional  (2-D) 
array  of  points  (xsk,  ysk;  k  =  1,2 , ,n).  Correspond¬ 
ing  to  illumination  of  the  Mh  grid  point  on  the  source 
plane,  a  charge-coupled  device  (CCD)  camera  records 
the  spatial  intensity  distribution,  Ik(xd,  yd),  on  the  exit 
surface  (or  detector  plane).  Every  pixel  of  the  CCD 
camera  thus  acts  as  a  detector  implementing  the  mul¬ 
tidetector  measurement  arrangement.  The  differ¬ 
ence,  AIk{xd,yd),  between  the  above-mentioned 
spatial  intensity  distribution,  Ik(xd,yd),  and  an  esti¬ 
mated  background  (say,  an  averaged  intensity  distri¬ 
bution  from  different  source  scanning  positions) 
provides  the  perturbation  in  the  spatial  intensity  dis¬ 
tribution  in  the  detector  plane  for  illumination  at  the 
&th  grid  point. 

The  localization  algorithm  is  based  on  the  premise 
that  each  object  (or  inhomogeneity)  within  the  turbid 
medium  alters  the  propagation  of  light  through  the 
medium.  Consequently,  the  spatial  distribution  of  the 
light  intensity  at  the  detector  plane  of  the  medium  is 
different  with  embedded  inhomogeneities  than  that 
without  them.  %he  influence  of  an  object  on  AIk(xd,  yd) 
involves  propagation  of  light  from  the  source  to  the 
object,  and  from  the  object  to  the  detector,  and  can  be 
described  in  terms  of  two  Green’s  functions  (propa¬ 
gators):  The  first  G( r,  rs)  describes  light  propagation 
from  the  source  at  rs  to  the  object  at  r,  and  the  second 


G(rd,  r)  describes  that  from  the  object  to  the  detector 
at  rd.  To  correlate  the  perturbations  in  the  light  in¬ 
tensity  distributions,  AIk(xd,  yd),  with  the  objects  em¬ 
bedded  in  the  turbid  medium,  we  assumed  that  these 
objects  illuminated  by  the  incident  wave  are  virtual 
sources  and  that  AIk(xd,yd)  are  taken  to  be  some 
weighted  mixtures  of  signals  arriving  from  these  vir¬ 
tual  sources  to  the  detector  plane.  ICA  assumes  these 
virtual  sources  to  be  independent,  and  based  on  that 
assumption  it  provides  the  independent  components. 
The  number  of  leading  independent  components  is 
the  same  as  the  number  of  the  embedded  objects.  The 
effective  contributions  of  independent  components  to 
the  light  intensity  distribution  on  the  source  and  de¬ 
tector  planes  are  proportional  to  the  projection  of  the 
Green’s  function,  G(r,  rs)  and  G(rd,  r),  on  the  source 
and  detector  planes,  respectively.  The  location  and 
characteristics  of  the  objects  are  obtained  from  fitting 
either  or  both  of  these  projections  to  those  of  the 
model  Green’s  function  in  the  background  medium. 

The  remainder  of  the  article  is  organized  as  follows. 
In  Section  2,  we  present  the  general  theoretical 
framework  for  OPTICA  and  then  discuss  the  specific 
case  of  a  turbid  medium  in  the  form  of  a  slab.  How¬ 
ever,  the  approach  can  be  adapted  to  any  arbitrary 
geometry.  Scattering  and  absorbing  objects  are  con¬ 
sidered  separately.  Section  3  presents  the  experimen¬ 
tal  methods,  materials,  and  parameters.  The  results 
are  presented  in  Section  4.  Finally,  the  implications 
of  these  results  and  the  scope  of  OPTICA  are  dis¬ 
cussed  in  Section  5. 


2.  Theoretical  Formalism 

In  the  linearized  scheme  of  inversion,  the  perturba¬ 
tion  of  the  detected  light  intensities  on  the  bound¬ 
aries  of  the  medium  (the  scattered  wave  field)  due  to 
absorptive  and  scattering  objects  (inhomogeneities) 
is  given  by3’13 

<kca(rd,  rs)  =  —  J  G(rd,  r)8p,a(r)cG(r,  rjd3r 

-J  d3r8D(r)cVrG(rd,  r)-VrG(r,  r„) 

(1) 

in  the  DA  (Ref.  31)  when  illuminated  by  a  unit  point 
source,  where  rs,  r  and  rrf  are  the  positions  of  the 
source,  the  inhomogeneity,  and  the  detector,  respec¬ 
tively;  8p.a  =  |xa  obj  -  p.a  and  8 D  =  Dobj  -  D  are  the 
differences  in  the  absorption  coefficient  and  the  dif¬ 
fusion  coefficient,  respectively,  between  the  inho¬ 
mogeneity  and  the  background;  c  is  the  speed  of 
light  in  the  medium;  and  G(r,  r')  is  the  Green’s  func¬ 
tion  describing  light  propagation  from  r'  to  r  inside 
the  background  turbid  medium  of  the  absorption  and 
diffusion  coefficients  |jl0  and  D,  respectively.  We  do 
not  explicitly  include  the  modulation  frequency  «  of 
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the  incident  wave  in  the  arguments  of  Eq.  (1)  for 
clarity.  The  following  formalism  can  be  applied 
to  continuous-wave,  frequency-domain,  and  time- 
resolved  measurements.  The  time-domain  measure¬ 
ment  is  first  Fourier  transformed  over  time  to  obtain 
data  over  many  different  frequencies.  Although  Eq. 
(1)  starts  with  DA,  it  should  be  emphasized  that  the 
formalism  is  not  limited  to  DA,  but  can  be  used  with 
other  models  of  light  propagation  in  turbid  media, 
such  as  the  cumulant  approximation,20-22'32  the 
random-walk  model,10-24  and  radiative  transfer17-33 
when  they  are  linearized. 

The  Green’s  function  G  for  a  slab  geometry  in  the 
diffusion  approximation  is  given  by 


G(r,  r')  =  G(p,  z,  z') 


'  4 ttD  kt- 


exp(- «r/)  exp(  — xr*  ) 


rk 


(2) 


rk±  =  [p2  +  (zTz'±2kd)2]1/2 


for  an  incident  amplitude-modulated  wave  of  modu¬ 
lation  frequency  <v,  where  k  -  0,  ±1,  ±2, . . . ,  p  = 
jjc  —  x')2  +  (y  -  y')2]1/2  is  the  distance  between  the  two 
points  r  =  (x,  y,  z)  and  r'  =  (x',y',  z')  projected  onto 
the  xy  plane,  k  =  [(p.a  -  iu>/c)/Df/2  chosen  to  have  a 
nonnegative  real  part,  and  the  extrapolated  bound¬ 
aries  of  the  slab  are  located  at  z  =  0  and  z  -  d 
=  L2  +  2 ze,  respectively,  where  Lz  is  the  physical 
thickness  of  the  slab  and  the  extrapolation  length  ze 
should  be  determined  from  the  boundary  condition  of 
the  slab.34-36  Equation  (2)  serves  as  the  model 
Green’s  function  for  the  uniform  background  medium 
of  a  slab  geometry.  The  modulation  frequency  to 
=  0  for  continuous-wave  light. 

In  practice,  the  projections  of  the  Green’s  function 
on  the  source  and  detector  planes  are  determined 
from  the  measured  perturbations  in  the  light  inten¬ 
sity  distribution  through  independent  component 
analysis.  The  comparison  with  the  Green’s  function 
computed  with  Eq.  (2)  is  then  used  to  locate  and 
characterize  the  inhomogeneities.  We  develop  the 
formalism  for  absorptive  and  scattering  inhomoge¬ 
neities  in  the  Subsections  2.A  and  2.B,  respectively. 

A.  Absorptive  Inhomogeneity 

The  assumption  that  absorptive  inhomogeneities  are 
localized  [that  is,  the jth  one  is  contained  in  volume  V, 
centered  at  r/l  J)]  enables  one  to  rewrite  the 

scattered  wave  field  in  Eq.  (1)  as 

j 

-4>sca(rrf,  rs)  =  2  G( rd,  r)qG(r-,  rs),  (3) 

7=1 

where  q:  =  8p,a(r7)cV,  is  the  absorption  strength  of  the 
jth  inhomogeneity.  The  scattered  wave  may  be  inter¬ 
preted  as  an  instantaneous  linear  mixture37 

x(rs)=As(rs).  (4) 


Here  s(rs)  =  [giG(r1;  rs), .  .  .  ,  qjG(rj,  r,)]T  represents 
the  J  virtual  sources;  i.e.,  the  J  inhomogeneities  illu¬ 
minated  by  the  incident  wave,  A  is  the  mixing  matrix 
given  by 


"G(rdl,  rj) 

G( rdl,  r2)  .  . 

•  G(rdl,  TjV 

A  = 

G( rd2, 

G( rd2,  r2)  . . 

■  G( rd2,  rj) 

_G(rdm,  rj) 

G(rdm,  r2)  .  . 

■  •  G(Yd'm,  Yj)_ 

whose  jth  column  (mixing  vector)  provides  the 
weight  factors  for  the  contributions  from  the  jth 
inhomogeneity  to  the  detectors,  and  x(rs)  = 
[-4>sca(rdl,  rs), . . .  ,  -c))sca(rdm,  rs)f  is  the  observed 
light  intensity  change.  The  superscript  T  denotes 
transposition.  The  incident  light  source  scans  a  total 
of  n  positions  rSl, . .  . ,  rSn  sequentially.  For  each 
source  position  r„  the  observation  is  made  over  m 
positions  rdl, . . . ,  rdm.  Each  set  of  such  measurement 
is  considered  data  at  one  temporal  sampling  point,  as 
used  in  the  conventional  instantaneous  linear  mix¬ 
ture  model.38  The  multisource,  multidetector  data  set 
x(rs)  thus  describes  signals  observed  in  m  channels 
( m  detectors)  from  J  virtual  sources  (or  J  inhomoge¬ 
neities)  simultaneously  over  n  discrete  temporal 
points  ( n  spatial  scanning  points).  One  absorptive 
inhomogeneity  is  represented  by  one  virtual  source 
qfiiYj,  rs).  The  virtual  source  qjG(Yj,  rs)  represents 
the  individual  inhomogeneity  illuminated  by  the  in¬ 
cident  wave  and  is  similar  to  the  concept  of  the  sec¬ 
ondary  source  in  Huygen’s  principle.39  The  role  of 
detectors  and  sources  can  be  interchanged  owing  to 
the  reciprocal  property  of  light  propagation. 

The  principal  assumption  of  this  formalism  is  that 
the  virtual  source  q,G(r,,  rs)  at  the  jth  inhomogeneity 
is  independent  of  the  virtual  sources  at  other  loca¬ 
tions.  Under  this  assumption,  ICA  can  be  used  with 
the  observations  from  the  light  source  scanned  at 
n  3>  J  positions  to  separate  out  both  virtual  sources 
s(r5)  and  the  mixing  matrix  A.26’37 

ICA  is  a  statistical  approach  to  separate  indepen¬ 
dent  sources  from  linear  instantaneous  or  convolu- 
tive  mixtures  of  independent  signals  without  relying 
on  any  specific  knowledge  of  the  sources  except  that 
they  are  independent.  The  sources  are  recovered  by  a 
minimization  of  a  measure  of  dependence,  such  as 
mutual  information,26-27  between  the  reconstructed 
sources.30-37  The  recovered  virtual  sources  and  mix¬ 
ing  vectors  from  ICA  are  unique  up  to  permutation 
and  scaling.30-37 

The  two  Green’s  functions  of  light  propagating 
from  the  source  to  the  inhomogeneity  and  from  the 
inhomogeneity  to  the  detector  are  retrieved  from  the 
separated  virtual  sources  s(rs)  and  the  mixing  matrix 
A.  The  jth  element  Sj(ys)  of  the  virtual  source  array 
and  the  jth  column  a,  (mixing  vector)  of  the  mixing 
matrix  A  provide  the  scaled  projections  of  the  Green’s 
function  on  the  source  and  detector  planes,  G(r7,  rs) 
and  G( rd,  rj),  respectively.  We  can  write 
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s/rs)  =  c tfiirj,  rs), 


aArd)  =  ^p(rd,  rj). 


medium,  and  Vj  is  the  volume  of  the  jth  scattering 
inhomogeneity,  the  scattered  wave  field  can  be  trans¬ 
formed  to 


where  a,-  and  (3,  are  scaling  constants  for  the  jth  in¬ 
homogeneity. 

Both  the  location  and  the  strength  of  the  jth  object 
can  be  computed  by  a  simple  fitting  procedure  by  use 
of  Eq.  (6).  We  adopted  a  least-square  fitting  procedure 
given  by 


j'  J' 

-4>sca(rd,  rs)  =  rd)q/gz(rj,  rs)  +  2  Pdj 

j= i 

X  cos  PdgJXj,  rd)q/psJ  cos  Osg,(r,,  rs) 

J' 

+  2  p dj  smOjg  irj,  rd)q/psj 
j=  i 

X  sin  %sgL{rj,  rs),  (11) 


min  |2  kj  1sj{rs)-G{rj,  rs)]2  +  2  IP,  ^X1^) 

-  G(rd,  r,)]2}.  (7) 

The  fitting  yields  the  location  r,  of  and  the  two  scaling 
constants  a j  and  for  the  jth  inhomogeneity,  whose 
absorption  strength  is  then  given  by  q}  =  afij. 

B.  Scattering  Inhomogeneity 

For  scattering  inhomogeneities,  under  the  assump¬ 
tion  that  the  inhomogeneities  are  localized  in  a  few 
regions,  the  same  analysis  can  be  carried  out  as  that 
for  absorptive  inhomogeneities.  The  only  modifica¬ 
tion  is  that  up  to  three  virtual  sources  may  appear  for 
one  scattering  inhomogeneity  corresponding  to  the 
x,  y,  z  components  in  the  dot  product  VrG(rrf,  r)  • 
VrG(r,  rs)  =  dxG( rd,  r)dxG(r,  rs)  +  dyG(rd,  r)dyG(r,  rs) 
+  32G( rd,  r)dzG(r,  rs)  in  Eq.  (1). 

Introducing  two  auxiliary  functions 


£Xr,  = 


(k  r*+  +  1) 


exp(-KrA') 


(ri+)3 


-  (kt*  +  1) 


exp(-Kr*  ) 


(ri  )' 


(8) 


1  ±2 


gz(r,  r ')  =  ^2  j(z  -z'  +  2kd)Urk+  +  1) 

exp(-«r*+) 


X 


(V)3 


X  (k rk  +  1) 


■  (z  +  z'  —  2  kd) 


exp(-K rk  ) 

(rk~)3 


(9) 


and  the  scattered  wave  due  to  scattering  inhomoge¬ 
neities  can  be  rewritten  as 


XscXri/,  rs)  =  —  J  dsr8ZXr)c{[(x  -  xd)(x  -  xs) 

+  (y-J'd)Cy-ys)]^i(r,  r d)g±(r,  rs) 

+  gz(r,  rd)gz(r,  rs)}.  (10) 

Denoting  the  scattering  inhomogeneities  as  q/ 
=  bD(rj)cVj ,  where  c  is  the  speed  of  light  in  the 


where  pdj  =  [( xd  -  xj)2  +  (yd  -  y/]1/2,  P v  =  [fe 
—  Xj)2  +  (ys  -  y,)2]1/2>  and  0d  and  0S  are  the  azimuthal 
angles  of  rd  —  rd  and  rs  -  Tj,  respectively.  This  scat¬ 
tered  wave  can  be  regarded  as  a  mixture  of  contribu¬ 
tions  from  (3 J')  virtual  sources: 

Qj'gJXj,  rs),  q/psJ  cos  $£  ,  (r,-,  rs), 

Q/Psj  sin  0.^  j  (r^  rs),  (12) 


with  mixing  vectors 

gJ-Xj,  rd),  pdj  cos  i)dgL(xP  rd),  pdj  sin  0dg±(r7-,  rd), 

(13) 

where  1  <j  <  J',  respectively.  There  are  in  general 
three  virtual  sources  of  specific  patterns  (one  cen- 
trosymmetric  and  two  dumbbell  shaped)  associated 
with  one  scattering  inhomogeneity,  whereas  only  one 
centrosymmetric  virtual  source  is  associated  with  one 
absorptive  inhomogeneity.  This  difference  may  be 
used  to  discriminate  absorptive  inhomogeneities 
from  scattering  inhomogeneities.  However,  for  scat¬ 
tering  inhomogeneities  deep  within  turbid  media, 
only  the  q/gz(rj,  rs)  virtual  source  remains  significant 
and  the  other  two  are  much  diminished.  In  such  a 
situation,  other  corroborative  evidence,  such  as  mul¬ 
tiwavelength  measurements,  are  required  to  deter¬ 
mine  the  nature  of  inhomogeneities.  Both  the 
location  and  the  strength  of  the  jth  scattering  object 
are  computed  by  fitting  the  retrieved  virtual  sources 
and  mixing  vectors  to  expressions  (12)  and  (13),  re¬ 
spectively. 

No  specific  light  propagation  model  is  assumed  in 
ICA.  The  only  assumption  is  that  virtual  sources  are 
mutually  independent.  The  number  of  inhomogene¬ 
ities  within  the  medium  is  determined  by  the  number 
of  the  independent  components  presented  in  the  mul¬ 
tisource,  multidetector  data  set.  The  analysis  of  re¬ 
trieved  independent  components  from  ICA  then 
localizes  and  characterizes  the  absorptive  and  scat¬ 
tering  inhomogeneities  inside  the  turbid  medium  in 
which  an  appropriate  model  of  the  light  propagator  is 
adopted.  When  the  noise  level  is  high  or  systematic 
errors  are  present,  or  both,  extra  independent  com¬ 
ponents  may  appear.  Only  the  leading  independent 
components  need  to  be  analyzed  to  detect  and  char¬ 
acterize  the  inhomogeneities  of  interest,  and  other 
components  can  be  discarded. 
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Fig.  1.  Schematic  diagram  of  specimen  1  comprising  an 
Intralipid-10%  suspension  in  water  with  two  long  cylindrical  ab¬ 
sorbing  objects  of  absorption  coefficient  0.23  mm-1. 


3.  Experimental  Methods  and  Materials 

Two  tissue-simulating  phantoms  with  absorption 
and  scattering  coefficients  within  the  reported  range 
of  values  for  healthy  human  breast  tissue  were  used 
in  the  study  reported  here.40 

The  first  specimen  (1),  shown  schematically  in  Fig. 
1,  was  a  250  mm  X  250  mm  X  50  mm  transparent 
plastic  container  filled  with  Intralipid-10%  suspen¬ 
sion  in  water.  The  concentration  of  Intralipid- 
10%  was  adjusted41  to  provide  a  transport  length  lt 
~  1  mm  and  an  absorption  coefficient  ga  = 
0.003  mm-1  at  785  nm,  emulating  those  of  human 
breast  tissue.  Two  cylindrical  glass  tubes  (outer  di¬ 
ameter,  8  mm;  inner  diameter,  6.98  mm;  and  length, 
250  mm)  were  filled  with  an  Intralipid-10%  suspen¬ 
sion  to  provide  the  same  scattering  coefficient,  but 
the  absorption  coefficient  was  changed  to  0.023  mm-1 
by  the  addition  of  absorbing  ink.  The  two  absorptive 
rods  are  placed  at  (x,  z)  =  (24,  29)  mm  and  (x,  z) 
=  (47,  33)  mm,  respectively,  with  the  fixes  of  cylin¬ 
drical  tubes  along  y. 

The  second  specimen  (2),  loaned  to  us  by  J.  C. 
Hebden  of  University  College  London  and  displayed 
schematically  in  Fig.  2,  was  a  166-mm-long,  82-mm- 
wide,  and  55-mm-thick  slab  made  of  materials  with  a 
reduced  scattering  coefficient  p,s'  ~  0.9  mnT1  (trans¬ 
port  length,  lt  ~  1.1  mm)  and  an  absorption  coeffi¬ 
cient  ~  0.006  mm-1.  The  slab  contained  four  5- 
mm-diameter,  5-mm-long  cylindrical  inhomogene¬ 
ities.  The  center  of  each  cylinder  was  located  in  the 
plane  halfway  between  the  front  and  the  back  faces 
of  the  slabs.  The  absorption  coefficient  of  each  cyl¬ 
inder  was  0.006  mm-1,  the  same  as  that  of  the  ma¬ 
terial  of  the  slab,  but  the  scattering  coefficients 
were  4,  2,  1.5,  and  1.1  times  greater.  The  first  and 
the  third  cylinders,  and  the  second  and  the  fourth 
cylinders,  are  on  two  horizontal  lines  approximately 


Fig.  2.  Schematic  diagram  of  specimen  2  obtained  from  Univer¬ 
sity  College  London.  It  is  a  solid  rectangular  block  embedded  with 
four  5-mm-diameter  and  5-mm-long  scattering  cylindrical  objects 
with  their  centers  on  the  central  plane.  The  absorption  and  scat¬ 
tering  characteristics  of  the  specimens  and  the  lateral  positions  of 
the  four  cylinders  are  described  in  the  text. 


22  mm  apart.  The  distance  between  neighboring 
cylinders  is  11  mm.  Further  details  about  the  slab 
may  be  obtained  from  an  article  published  by  Hall 
et  al42 

The  experimental  arrangement  used  for  imaging 
of  these  two  specimens  (1  and  2)  is  shown  schemat¬ 
ically  in  Fig.  3.  For  cw  measurements,  a  200-pm 
fiber  delivered  a  beam  of  784-nm  light  from  a  diode 
laser  (Ocean  Optics  R-2000)  to  illuminate  the  input 
surface  (or  source  plane)  of  the  specimen.  A  cooled 
CCD  camera  set  at  an  acquisition  time  of  150  ms 
recorded  2-D  intensity  patterns  of  the  light  trans¬ 
mitted  through  the  opposite  side  of  the  slab  speci¬ 
men.  For  time-resolved  measurements,  we  used  a 
1-mm-diameter  collimated  beam  of  784-nm,  150-fs, 
1-kHz  repetition  rate  light  pulses  from  a  Ti:sap- 
phire  laser  and  amplifier  system43  for  sample  illu- 


Source 

plane 


Detector 

plane 


Fig.  3.  Schematic  diagram  of  the  experimental  arrangement  for 
imaging  objects  embedded  in  a  turbid  medium.  Inset  shows  the  2-D 
array  in  the  input  plane  that  is  scanned  across  the  incident  laser 
beam. 
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mination.  An  ultrafast  gated  intensified  camera 
system  (UGICS),  which  provided  a  FWHM  gate 
width  variable  from  80  ps  to  6  ns,  recorded  the  2-D 
intensity  patterns  of  the  light  transmitted  through 
the  opposite  side  of  the  slab.  Computer-controlled 
xy  translation  stages  scanned  the  specimens  in  an 
array  of  points  in  the  xy  plane,  as  displayed  in  Fig.  3. 
For  the  long  cylindrical  tubes  in  specimen  1,  a  line 
scan  of  16  points  with  a  step  size  of  2.5  mm  along  the 
x  axis  was  enough  to  obtain  the  (x,  z)  locations  of  the 
absorbing  cylinders.  An  array  of  20  X  18  points  with 
a  step  size  of  2.5  mm  across  the  lateral  positions  of 
the  four  inhomogeneities  in  specimen  2  was 
scanned  to  obtain  their  locations. 


4.  Results 

Temporal  profiles  of  the  transmitted  pulses  were 
measured  by  use  of  the  UGICS  in  the  scan  mode  with 
an  80-ps  gate  width.  The  average  optical  properties  of 
the  turbid  medium  were  estimated  by  fitting  the  tem¬ 
poral  profiles  to  the  DA  of  the  RTE  for  a  slab  geometry. 

ICA  of  the  perturbations  in  the  spatial  intensity 
distributions  provided  the  corresponding  indepen¬ 
dent  intensity  distributions  on  the  source  and  de¬ 
tector  planes.  ICA-generated  independent  intensity 
distributions  on  the  source  and  detector  planes  are 
shown  in  the  first  and  second  rows,  respectively,  of 
Fig.  4  for  the  two  absorbing  cylinders  in  specimen  1. 
The  locations  of  the  absorbing  cylinders  are  ob¬ 
tained  from  fitting  these  independent  component 
intensity  distributions  to  those  of  the  DA  in  a  slab 
[Eq.  (2)] .  The  first  cylinder  is  found  at  x  =  24  mm, 
29  mm  away  from  the  source  plane  and  21  mm  away 
from  the  detector  plane.  The  second  cylinder  is  found 
at  x  =  47  mm,  33  mm  away  from  the  source  plane 
and  17  mm  away  from  the  detector  plane.  The  (x,  z) 
coordinates  of  both  the  cylinders  agree  to  within 
0.5  mm  of  their  known  locations.  The  absorption 
strengths  of  the  two  rods  are  estimated  by  use  of  a 
least-square  fitting  procedure  [Eq.  (7)].  The  resolved 
absorption  strengths  are  qx  =  0.152  mm2/ps  and  q2 
=  0.132  mm2/ps,  respectively,  for  the  left  and  right 
rods.  The  values  are  88%  and  76%,  respectively,  of 
the  true  value  of  q  =  0.173  mm2/ps. 

The  independent  intensity  distributions  at  the  de¬ 
tector  plane  corresponding  to  the  four  scattering  inho¬ 
mogeneities  in  specimen  2  are  displayed  in  Figs.  5(a)- 
5(d).  These  independent  components  are  then  used  to 
obtain  the  projections  of  the  inhomogeneity  detector 
Green’s  function,  G(rd,  r ,),j  =  1,  2,  3,  4,  on  the  detec¬ 
tor  plane  for  the  four  small  cylindrical  scattering  in¬ 
homogeneities  embedded  in  specimen  2.  The 
locations  of  the  inhomogeneities  are  determined  by 
fitting  the  projections  to  those  of  the  model  Green’s 
function.  The  locations  of  all  four  inhomogeneities 
were  obtained.  Even  the  weakest  scatterer,  with  a 
scattering  coefficient  just  1.1  times  the  background 
and  hence  considered  to  be  rather  unlikely  to  be 
found,42  was  detected.  Positions  along  the  z  axis 
(depth)  of  the  cylinders  were  found  to  be  at  28.1,  27.9, 
27.1,  and  32.6  mm.  Except  for  the  last  cylinder,  the 


Fig.  4.  Normalized  independent  spatial  intensity  distributions  as 
a  function  of  the  lateral  position  x  at  the  input  (or  source)  plane 
(first  row)  and  the  exit  (or  detector)  plane  (second  row)  generated 
by  ICA  for  the  two  absorbing  cylinders  in  specimen  1.  The  hori¬ 
zontal  profile  of  the  intensity  distributions  on  the  source  plane 
(diamond)  and  on  the  detector  plane  (circle)  are  displayed  in  the 
third  row.  Solid  curves  show  the  respective  Green’s-function  fit 
used  for  obtaining  the  locations  of  the  objects. 


depth  of  the  cylinders  agree  well  with  their  known 
center  positions  of  27.5  mm.  The  lateral  positions  are 
determined  to  be  (62,  63),  (48,  33),  (33,  62),  and 
(18,  33)  mm  for  the  four  scattering  cylinders  (see  Ta¬ 
ble  1).  The  strongest  and  the  third-strongest  scatter- 
ers  are  on  the  same  horizontal  line  y  —  62  mm, 
whereas  the  second-strongest  and  the  weakest  scat¬ 
tered  are  on  the  horizontal  line  y  —  33  mm  with  a 
spacing  of  29  mm.  The  four  scatterers  are  separated 
by  equal  spacing,  —14  mm  in  the  horizontal  direc¬ 
tion.  The  lateral  positions  agree  well  with  the  known 
(x,  y)  coordinate  values.  The  uncertainties  in  location 
and  separation  are  not  greater  than  3  mm  except  for 
the  weakest  target. 

5.  Discussion 

The  OPTICA  presented  in  this  article  introduces  the 
information  theory  technique  of  ICA  to  the  problem  of 
optical  tomographic  imaging  of  objects  in  turbid  me¬ 
dia.  It  is  shown  to  provide  object  locations  accurate  to 
— 1  mm  in  human-breast-like  turbid  media.  It  uses 
multiple-source  (realized  in  this  case  through  scan¬ 
ning  of  the  sample  in  the  xy  plane  across  the  incident 
beam  propagating  in  the  z  direction)  illumination  and 
a  multiple-detector  (each  pixel  on  the  CCD  may  be 
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Fig.  5.  Independent  spatial  intensity  distributions  at  the  exit  (or 
detector)  plane  generated  by  ICA  corresponding  to  objects  with 
scattering  coefficients:  (a)  4  times,  (b)  2  times,  (c)  1.5  times,  and  (d) 
1.1  times  that  of  the  material  of  the  slab  in  specimen  2.  Horizontal 
profiles  of  the  intensity  distributions  in  (aHd)  are  shown  by  circles 
in  (e)-(h),  respectively,  with  solid  curves  representing  the  Green’s- 
function  fit  used  for  extracting  object  locations. 


viewed  as  a  detector)  data-acquisition  scheme.  The 
resulting  spatial  diversity  and  multiple  angular  ob¬ 
servations  provide  robust  data  for  extracting  3-D  lo¬ 
cation  information  about  the  embedded  objects 
(inhomogeneities)  in  the  medium.  A  salient  feature  of 
OPTICA  is  that  ICA  provides  the  independent  com¬ 
ponents  due  to  the  inhomogeneities  with  minimal 


Table  1.  Comparison  of  Known  and  OPTICA-Determined  Positions  of 
the  Four  Targets0 


Target 

Number 

Target 

Strength 

Known  Position 
(x,  y,  z)  (mm) 

OPTICA-Estimated 
Position 
(x,  y,  z)  (mm) 

1 

4 

(60,  60,  27.5) 

(62,  63,  28.1) 

2 

2 

(47,  30,  27.5) 

(48,  33,  27.9) 

3 

1.5 

(33,  60,  27.5) 

(33,  62,  27.1) 

4 

1.1 

(20,  30,  27.5) 

(18,  33,  32.6) 

"Target  strength  is  the  ratio  of  the  scattering  coefficients  of  the 
target  to  that  of  the  surrounding  medium.  The  errors  in  location 
are  not  greater  than  3  mm. 


processing  of  the  data  and  does  not  have  to  resort  to 
any  specific  light  propagation  model  for  obtaining 
this  information.  Specific  light  propagation  models 
are  needed  only  in  the  later  stage  to  determine  the 
location  by  curve  fitting  of  the  Green’s  functions. 
OPTICA  is  not  model  specific;  any  appropriate 
model  for  light  propagation,  including  the  DA  and 
the  cumulant  solutions  of  the  RTE,  may  be  used. 
OPTICA  can  be  used  with  contrast  agents  such  as 
fluorescence-based  optical  tomography  as  well. 

Although  we  used  the  slab  geometry  in  the  study 
reported  in  this  article,  the  approach  does  not  depend 
on  any  specific  geometry.  It  may  be  used  for  other 
geometries  or  even  an  arbitrary-shaped  boundary.  The 
approach  is  fast  and  is  expected  to  be  amenable  to 
near-real-time  detection  and  localization  of  objects  in  a 
turbid  medium,  which  is  a  key  consideration  for  in  vivo 
medical  imaging.  The  approach  is  remarkably  sensi¬ 
tive,  considering  that  it  could  discern  all  four  cylinders 
in  specimen  2.  The  approach  successfully  detected 
even  the  lowest-contrast  inhomogeneity  of  the  four 
that  had  a  reduced  scattering  coefficient  only  10% 
higher  than  the  surrounding  medium  and  was  consid¬ 
ered  improbable  to  be  detected.42  OPTICA  obtains  lo¬ 
cations  of  the  objects  by  fitting  either  or  both  of  the 
Green’s  functions  G(r,  rs)  and  G(rd,  r),  and  is  suited 
for  physically  small  inhomogeneities.  Given  its  ca¬ 
pability  of  identifying  low-contrast  small  objects, 
the  approach  is  expected  to  be  useful  for  detecting 
tumors  at  their  early  stages  of  development,  a  cov¬ 
eted  goal  in  medical  imaging. 

As  demonstrated  with  specimen  1  and  specimen  2, 
the  approach  could  locate  both  absorptive  and  scat¬ 
tering  objects.  When  both  absorptive  and  scattering 
objects  are  present  in  the  same  turbid  specimen,  OP¬ 
TICA  can  locate  them,  but  their  identification  as  ab¬ 
sorbing  or  scattering  entities  becomes  a  more 
challenging  task.  As  discussed  in  connection  with  ex¬ 
pressions  (12)  and  (13),  each  scattering  inhomogene¬ 
ity  is  expected  to  be  represented  by  three  virtual 
sources,  yielding  three  pairs  of  effective  intensity  dis¬ 
tributions  each  on  the  detector  and  source  planes. 
When  the  background  scattering  is  not  severe  and 
the  S/N  ratio  is  high,  contributions  from  all  three 
virtual  sources  may  be  distinguished  and  the  corre¬ 
sponding  object  may  be  identified  as  a  scattering  en¬ 
tity.  Our  simulation  results  support  this  assertion. 
For  highly  scattering  conditions  with  lower  S/N  ra¬ 
tios,  contributions  from  the  two  dumbbell-shaped  vir¬ 
tual  sources  may  not  be  discerned,  and  corroborative 
information  obtained  by  other  means,  such  as  mea¬ 
surements  using  light  of  different  wavelengths,  are 
required  for  identification  of  the  object  as  an  absorp¬ 
tive  or  scattering  entity.  Multiwavelength  spectro¬ 
scopic  imaging  measurements  have  the  potential  to 
provide  diagnostic  information,  such  as  whether  a 
tumor  is  malignant  or  benign. 

In  summary,  OPTICA  has  the  potential  to  emerge 
as  a  new  versatile  tool  for  locating  targets  in  turbid 
media,  particularly  in  diagnostic  medical  imaging 
and  underwater  imaging. 
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Optical  imaging  of  turbid  media  using  independent 
component  analysis:  Theory  and  Simulation 
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Abstract 

A  new  imaging  approach  for  three-dimensional  localization  and  characterization  of  objects 
in  a  turbid  medium  using  independent  component  analysis  (ICA)  from  information  theory  is 
developed  and  demonstrated  using  simulated  data.  This  approach  uses  a  multi-source  and 
multi-detector  signal  acquisition  scheme.  Independent  component  analysis  of  the  perturba¬ 
tions  in  the  spatial  intensity  distribution  measured  on  the  medium  boundary  sorts  out  the 
embedded  objects.  The  locations  and  optical  characteristics  of  the  embedded  objects  are  ob¬ 
tained  from  a  Green’s  function  analysis  based  on  any  appropriate  model  for  light  propagation 
in  the  background  medium.  This  approach  is  shown  to  locate  and  characterize  absorptive  and 
scattering  inhomogeneities  within  highly  scattering  medium  to  a  high  degree  of  accuracy.  In 
particular,  we  show  this  approach  can  discriminate  between  absorptive  and  scattering  inhomo¬ 
geneities,  and  can  locate  and  characterize  complex  inhomogeneities  which  is  both  absorptive 
and  scattering.  The  influence  of  noise  and  uncertainty  in  background  absorption  or  scattering 
on  the  performance  of  this  approach  is  investigated. 

Keywords:  image  processing,  image  reconstruction,  medical  imaging,  inverse  problems,  absorption, 
scattering,  diffusion,  radiative  transfer 
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Optical  probing  of  the  interior  of  multiply  scattering  colloidal  suspensions  and  biological  materials 
has  attracted  much  attention  over  the  last  decade.  In  particular,  biomedical  optical  tomography 
and  spectroscopy  which  has  the  potential  to  provide  functional  information  about  brain  activities 
and  diagnostic  information  about  tumors  in  breast  and  prostate  are  being  actively  pursued[l,  2,  3, 
4,  5,  6,  7,  8,  9,  10,  11,  12,  13,  14,  15,  16,  17].  Simultaneous  developments  in  experimental  apparatus 
and  techniques  for  object  interrogation  and  signal  acquisition, [4,  5,  2,  18,  19]  analytical  models 
for  light  propagation,  [10,  20,  21,  22]  and  computer  algorithms  for  image  reconstruction [8,  6]  hold 
promise  for  realization  of  these  potentials  of  optical  tomography. 

Multiple  scattering  of  light  in  thick  turbid  media  precludes  direct  imaging  of  embedded  targets. 
One  typically  uses  an  inverse  image  reconstruction  (HR)  [8,  6]  approach  to  reconstruct  a  map  of  the 
optical  properties,  such  as,  absorption  coefficient  (fia)  and  scattering  coefficient  (/rs),  of  the  medium 
by  matching  the  measured  light  intensity  distribution  on  the  boundary  of  the  turbid  medium  to  that 
calculated  by  a  forward  model  for  the  propagation  of  light  in  that  medium.  The  commonly  used 
forward  models  include  the  radiative  transfer  equation  (RTE),[23,  17]  the  diffusion  approximation 
(DA)  of  RTE,[8,  6]  and  random  walk  of  photons.  [24,  25] 

The  inversion  problem  is  ill-posed  and  needs  to  be  regularized  to  stabilize  the  inversion  at  a 
cost  of  reduced  resolution.  [26,  8]  Both  iterative  reconstruction  and  noniterative  linearized  inversion 
approaches  have  been  used  to  solve  the  inversion  problem  in  optical  tomography,  which  is  weakly 
nonlinear,  with  limited  success.  Reconstruction  of  images  with  adequate  spatial  resolution  and 
accurate  localization  and  characterization  of  the  inhomogeneities  remain  a  formidable  task.  Time 
required  for  data  acquisition  and  image  reconstruction  is  another  important  consideration. 

In  this  article,  we  present  a  novel  algorithm  based  on  the  independent  component  analysis 
(ICA)[27,  28]  from  information  theory  to  locate  absorptive  and  scattering  inhomogeneities  embed¬ 
ded  in  a  thick  turbid  medium  and  demonstrate  the  efficacy  using  simulated  data.  ICA  has  been 
successfully  applied  in  various  other  application  such  as  electroencephalogram  and  nuclear  mag¬ 
netic  resonance  spectroscopy [29,  30,  28,  31]  We  refer  to  this  information  theory-inspired  approach 
as  optical  imaging  using  independent  component  analysis,  abbreviated  as,  OPTICA.  The  novelty 
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of  OPTICA  over  other  ICA  applications  is  that  OPTICA  associates  directly  the  independent  com¬ 
ponents  to  the  Green’s  functions  responsible  for  light  propagation  in  the  turbid  medium  from  the 
inhomogeneities  to  the  source  and  the  detector,  and  therefore  the  retrieved  independent  components 
can  be  used  to  locate  and  characterize  the  inhomogeneities. 

OPTICA  uses  a  multi-source  illumination  and  multi-detector  signal  acquisition  scheme  providing 
a  variety  of  spatial  and  angular  views  essential  for  three-dimensional  (3D)  object  localization.  Each 
object  (or,  inhomogeneity)  within  the  turbid  medium  alters  the  propagation  of  light  through  the 
medium.  The  influence  of  an  object  on  the  spatial  distribution  of  the  light  intensity  at  the  detector 
plane  involves  propagation  of  light  from  the  source  to  the  object,  and  from  the  object  to  the  detector, 
and  can  be  described  in  terms  of  two  Green’s  functions  (propagators)  describing  light  propagation 
from  source  to  the  object  and  that  from  the  object  to  the  detector,  respectively. 

The  absorptive  or  scattering  inhomogeneities  illuminated  by  the  incident  wave  are  assumed  to 
be  virtual  sources,  and  the  perturbation  of  the  spatial  distribution  of  the  light  intensity  on  the 
medium  boundary  is  taken  to  be  some  weighted  mixture  of  signals  arriving  from  these  virtual 
sources.  These  virtual  sources  are  statistically  independent  and  can  be  recovered  by  independent 
component  analysis  of  the  recorded  data  set.  The  number  of  leading  independent  components  is 
same  as  the  number  of  embedded  objects.  The  location  and  characterization  of  inhomogeneities 
are  obtained  from  the  analysis  of  the  retrieved  virtual  sources  using  an  appropriate  model  of  light 
propagation  in  the  background  medium. 

The  remainder  of  this  article  is  organized  as  follows.  In  Section  2,  we  present  a  brief  introduction 
to  independent  component  analysis  and  review  the  general  theoretical  framework  for  OPTICA. 
Section  3  presents  the  results  from  simulations  for  different  configurations.  Implications  of  these 
results  are  discussed  in  Section  4. 
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2  Theoretical  Formalism 


2.1  Independent  component  analysis 

Blind  source  separation  is  a  class  of  problem  of  general  interest  which  consists  of  recovering  un¬ 
observed  signals  or  virtual  sources  from  several  observed  mixtures.  Typically  the  observations  are 
the  output  of  a  set  of  sensors,  where  each  sensor  receives  a  different  combination  of  the  source 
signals.  Prior  knowledge  about  the  mixture  in  such  problems  is  usually  not  available.  The  lack  of 
prior  knowledge  is  compensated  by  a  statistically  strong  but  often  physically  plausible  assumption 
of  independence  between  the  source  signals.  Over  the  last  decade,  independent  component  analysis 
(ICA)  has  been  proposed  as  a  solution  to  the  blind  source  separation  problem  and  emerged  as  a 
new  paradigm  in  signal  processing  and  data  analysis. [27,  28,  32,  31] 

The  simplest  ICA  model,  an  instantaneous  linear  mixture  model  [32],  assumes  the  existence 
of  n  independent  signals  sl(t)  ( i  =  1,2, ...,n)  and  the  observation  of  at  least  as  many  mixtures 
Xi(t )  (i  —  1, 2, . . . ,  m)  by  m  >  n  sensors,  these  mixtures  being  linear  and  instantaneous,  i.e,  Xi(t)  = 
Y^j=i  aijSj(t)  for  each  i  at  a  sequence  of  time  t.  In  a  matrix  notation, 

x(/)  =  As{t),  A  G  Rmx"  (1) 

where  A  is  the  mixing  matrix.  The  jth  column  of  A  gives  the  mixing  vector  for  the  jth  virtual  source. 
Independent  component  analysis  can  be  formulated  as  the  computation  of  an  n  x  m  separating 
matrix  B  whose  output 

y(i)  =  Bx(i)  =  Cs(t),  B  G  R”xm,  C  =  BA  G  RnXn  (2) 

is  an  estimate  of  the  vector  s (t)  of  the  source  signals. 

The  basic  principle  of  ICA  can  be  understood  in  the  following  way.  The  central  limit  theorem 
in  probability  theory  tells  that  the  distribution  of  independent  random  variables  tends  toward 
a  Gaussian  distribution  under  certain  conditions.  Thus  a  sum  of  multiple  independent  random 
variables  usually  has  a  distribution  that  is  closer  to  Gaussian  than  any  of  the  original  random 
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variables,  y,(t)  =  Y jCijSj(t )  in  Eq.  2,  as  a  summation  of  independent  random  variables  Sj(t),  is 
usually  more  Gaussian  than  Sj(t)  while  y^t )  becomes  least  Gaussian  when  it  in  fact  equals  one  of 
the  Sj(t).  This  heuristic  argument  shows  that  independent  component  analysis  can  be  intuitively 
regarded  as  a  statistical  approach  to  find  the  separating  matrix  B  such  that  yt(i)  is  least  Gaussian. 
This  can  be  achieved  by  maximizing  some  measure  of  nongaussianity  such  as,  maximizing  kurtosis 
(the  fourth  order  cumulant),  of  yi{t).[ 33,  32] 

Independent  component  analysis  separates  independent  sources  from  linear  instantaneous  or 
convolutive  mixtures  of  independent  signals  without  relying  on  any  specific  knowledge  of  the  sources 
except  that  they  are  independent.  The  sources  are  recovered  by  a  maximization  of  a  measure  of 
independence  (or,  a  minimization  of  a  measure  of  dependence),  such  as  nongaussianity  and  mutual 
information  between  the  reconstructed  sources. [32,  31]  The  recovered  virtual  sources  and  mixing 
vectors  from  ICA  are  unique  up  to  permutation  and  scaling.  [32,  31] 

2.2  Optical  imaging  using  independent  component  analysis 

The  classical  approach  to  propagation  of  multiply  scattered  light  in  turbid  media,  which  assumes 
that  phases  are  uncorrelated  on  scales  larger  than  the  scattering  mean  free  path  ls ,  leads  to  the 
radiative  transport  equation  (RTE)  in  which  any  interference  effects  are  neglected.  [34]  RTE  does 
not  admit  closed-form  analytical  solutions  in  bounded  regions  and  its  numerical  solution  is  com¬ 
putational  expensive.  The  commonly  used  forward  models  in  optical  imaging  of  highly  scattering 
media  is  the  diffusion  approximation  (DA)  to  RTE.  [8,  6] 

The  approach  OPTICA  can  be  applied  to  different  models  of  light  propagation  in  turbid  media, 
such  as,  the  diffusion  approximation,^,  6]  the  cumulant  approximation[20,  35,  22],  the  random  walk 
model,  [24,  10]  and  radiative  transfer [34,  17]  when  they  are  linearized.  The  diffusion  approximation 
is  valid  when  the  inhomogeneities  are  deep  within  a  highly  scattering  medium.  We  only  discuss  the 
formalism  of  OPTICA  in  the  diffusion  approximation  in  this  article. 

In  the  diffusion  approximation,  the  perturbation  of  the  detected  light  intensities  on  the  bound¬ 
aries  of  the  medium,  the  scattered  wave  field,  due  to  absorptive  and  scattering  objects  (inhomo- 
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geneities)  can  be  written  as:  [3,  13] 


$Si 


x(rd,  r s)  =  ~  J  G(rd,  r)Spa(r)cG(r,  rs)dr  -  J  8D(r)cVTG(rd,  r)  •  VrG( r,  rs)dr 


(3) 


to  the  first  order  of  Born  approximation [36]  when  illuminated  by  a  point  source  of  unit  power, 
where  rs,  r  and  rd  are  the  positions  of  the  source,  the  inhomogeneity  and  the  detector  respectively, 
5p,„  =  pa,obj  —  na  and  8D  =  Z)ohj  —  D  are  the  differences  in  absorption  coefficient  and  diffusion 
coefficient,  respectively,  between  the  inhomogeneity  and  the  background,  c  is  the  speed  of  light  in 
the  medium,  and  G( r,  r')  is  the  Green’s  function  describing  light  propagation  from  r'  to  r  inside 
the  background  turbid  medium  of  absorption  and  diffusion  coefficients  pa  and  D. 

Eq.  (3)  is  written  in  the  frequency  domain  and  does  not  explicitly  show  the  modulation  fre¬ 
quency  oj  of  the  incident  wave  for  clarity.  The  following  formalism  applies  to  continuous  wave, 
frequency-domain  and  time-resolved  measurements.  The  time  domain  measurement  is  first  Fourier 
transformed  over  time  to  obtain  data  over  many  different  frequencies. 

The  Green’s  function  G  for  a  slab  geometry  in  DA  is  given  by 


G(r,r'; 


G{p,  z,  z')  = 


47 xD 


E 

k=-oo 


exp(-Kr+)  exp  (-«xfc)' 


(4) 


=  \] p2  +  (zTz'±  2  kd) 


for  an  incident  amplitude-modulated  wave  of  modulation  frequency  cu,  where  k  =  0,  ±1,±2, ..., 
p  =  y/\x  —  x ')2  +  (y  —  y')2  is  the  distance  between  the  two  points  r  =(x,y,z )  and  r'  =  (x',y',z') 
projected  onto  the  xy  plane,  n  —  y7 {pa  —  iuj/c)/ D  chosen  to  have  a  nonnegative  real  part,  and  the 
extrapolated  boundaries  of  the  slab  are  located  at  z  =  0  and  z  =  d  —  Lz  +  2 ze,  respectively,  where 
Lz  is  the  physical  thickness  of  the  slab  and  the  extrapolation  length  ze  should  be  determined  from 
the  boundary  condition  of  the  slab.  [37,  38,  39]  Greens’  functions  in  Eq.  (3)  for  other  geometries  can 
be  obtained  either  analytically  or  numerically.  [40,  41] 
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2.2.1  Absorptive  inhomogeneities 

We  first  consider  absorptive  inhomogeneities.  Under  the  assumption  that  absorptive  inhomo¬ 
geneities  are  localized,  that  is,  the  jth  one  is  contained  in  volume  Vj  centered  at  r,  (1  <  j  <  J), 
the  scattered  wave  field  in  Eq.  (3)  can  be  rewritten  as 

j 

-^sca(rd,rs)  =  Y^G{rd,rj)qjG(rj,rs)  (5) 

3= i 

where  qj  =  5jia(rj)cVj  is  the  absorption  strength  of  the  jth  inhomogeneity.  The  scattered  wave  is 
in  a  form  of  in  an  instantaneous  linear  mixture  (1).  One  absorptive  inhomogeneity  is  represented 
by  one  virtual  source  qjG(rj,rs )  with  a  mixing  vector  G(rd,  r^-). 

As  the  virtual  source  qjG(rj,  rs)  at  the  jth  inhomogeneity  is  independent  of  the  virtual  sources 
at  other  locations,  independent  component  analysis  can  be  used  with  the  observations  obtained 
for  the  light  source  at  n  »  J  different  positions  to  separate  out  both  virtual  sources  Sj( rs)  and 
the  mixing  vectors  aj(rd).  [27,  32]  The  jth  virtual  source  Sj(rs)  and  the  jth  mixing  vector  aj( rd) 
provide  the  scaled  projections  of  the  Green’s  function  on  the  source  and  detector  planes,  G(rj,rs) 
and  G(rd,rj),  respectively.  We  can  write 


Sj(  rs)  =  otjGi  rj,rs),  (6) 

aj{rd)  =  (3jG(rd,rj) 


where  otj  and  (5j  are  scaling  constants  for  the  jth  inhomogeneity. 

Both  the  location  and  strength  of  the  jth  object  can  be  computed  by  a  simple  fitting  procedure 
using  Eq.  (6).  We  adopted  a  least  square  fitting  procedure  given  by: 


min«  |  zL  ia3  lsJ'(rs)  _  G(rJ>rs)]2  +  Y,  [Pj  ‘fljW  -  G(r,j,  rj)] 2 1 

Tj  ,Oj  ,Pj  I 

t  r.i  rd  ) 


(7) 


The  fitting  yields  the  location  ij  of  and  the  two  scaling  constants  otj  and  /5j-  for  the  jth  inhomogeneity 
whose  absorption  strength  is  then  given  by  q3  —  otjfij. 


8 


2.2.2  Scattering  inhomogeneities 


For  scattering  inhomogeneities,  a  similar  analysis  shows  the  scattered  wave  can  be  written  as: 

-4>sca(rd,rs)  =  Y^9z(rj,rdWjgz(rj,ra)  (8) 

3=1 

r 

+  ^  Pdj  cos  0d9±(ij,  Yd)q'jPsj  COS  dsg±(rj,  r.s) 
i=i 
J' 

+  ^  Pdj  sin^i(r3)  rd)<?'psj  sin0s0X(r,-,  rs) 
i=i 


where  =  £Z)(rj)cVj  is  the  diffusion  strength  of  the  jth  scattering  inhomogeneity  of  volume  V- 
0  =  1,2,...,  J'),  Pdj  =  \/(xd  ~  x.j)2  +  (yd  -  r/j)2,  psj  =  -  x,)2  +  (?ps  -  y,)2>  #d  and  6»s  are  the 

azimuthal  angles  of  rd  —  and  rs  —  respectively,  and  the  two  auxiliary  functions  are  given  by 


5i.(rV)  = 


47tD 


+oo 

E 

k=— oo 


(Kr+  +  1) 


exp(- 


■  KT 


(rkY 


d_(sr- + 


-nr. 


—  \  3 


( rk ) 


(9) 


and 


1  +00  f 

9z{r,  ^  =  4 ^  j  “  z'  +  2fcc0(Krfc  +  1) 

k=— oo  V 


exp(-Kr^) 

K)3 


-(z  +  /-2fcd)(Krfc-  +  l)eXp(  ^T) 

(rJ 


(10) 


The  scattered  wave  from  one  scattering  inhomogeneity  is  thus  a  mixture  of  contributions  from 
(3  J')  virtual  sources: 


Qj9z(r3,rs),  q'jpsj  cos0spx(rj,  rs),  g'psjsin0sp±(rj,rs), 


(11) 


with  mixing  vectors 


,9z(rj)  r(;),  PdjCosedg±(vj,Yd),  prfj  sin  0dpx  (r, ,  rd) 


(12) 
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where  1  <  j  <  ,/',  respectively.  Both  the  location  and  strength  of  the  jth  scattering  object  are 
computed  by  fitting  the  retrieved  virtual  sources  and  mixing  vectors  to  Eq.  (11)  and  Eq.  (12)  using 
a  least  square  procedure,  respectively. 

There  are  in  general  three  virtual  sources  of  specific  patterns  (one  centrosymmetric  and  two 
dumbbell-shaped)  associated  with  one  scattering  inhomogeneity,  whereas  only  one  centrosymmetric 
virtual  source  is  associated  with  one  absorptive  inhomogeneity.  This  difference  may  serve  as  the 
basis  to  discriminate  absorptive  and  scattering  inhomogeneities. 

The  only  assumption  made  in  OPTICA  is  that  virtual  sources  are  mutually  independent.  The 
number  of  inhomogeneities  within  the  medium  is  determined  by  the  number  of  the  independent 
components  presented  in  the  multi-source  multi-detector  data  set.  No  specific  light  propagation 
model  is  assumed  in  this  step.  The  analysis  of  retrieved  independent  components  from  ICA  then 
localizes  and  characterizes  the  absorptive  and  scattering  inhomogeneities  inside  the  turbid  medium 
using  an  appropriate  model  of  the  light  propagator.  Extra  independent  components  may  appear 
depending  on  the  level  of  noise  in  the  data.  These  components  can  be  discarded  and  only  the 
leading  independent  components  need  to  be  analyzed  to  detect  and  characterize  inhomogeneities  of 
interest. 

3  Results 

Simulations  were  performed  for  a  slab  of  thickness  50mm  shown  schematically  in  Fig.  1.  The 
absorption  and  diffusion  coefficients  of  the  uniform  slab  is  a  =  l/300mm-1and  D  =  l/3mm 
respectively,  close  to  that  of  human  breast  tissue.  [42]  The  incident  CW  beam  scans  a  set  of  21  X  21 
grid  points  covering  an  area  of  90  x  90mm2.  The  spacing  between  two  consecutive  grid  points  is 
4.5mm.  The  light  intensity  on  the  other  side  of  the  slab  is  recorded  by  a  CCD  camera  on  42  x  42 
grid  points  covering  the  same  area. 


[Figure  1  about  here.] 

In  the  simulations  presented  in  the  following  subsections,  we  fix  the  ratio  of  strength  of  absorption 
to  that  of  diffusion  to  be  0.01,  which  produce  perturbations  of  comparable  magnitude  on  the  light 
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intensities  measured  on  the  detector  plane  from  the  absorption  and  scattering  inhomogeneities.  As 
the  scattered  wave  is  linear  with  respect  to  the  absorption  and  diffusion  strengths,  we  also  set  the 
strength  of  either  absorption  or  diffusion  to  be  unity  in  simulations  for  convenience. 

3.1  Absorptive  Inhomogeneities 

Two  absorptive  inhomogeneities,  each  of  a  unity  absorption  strength,  are  placed  at  positions 
(50, 60, 20)mm  and  (30, 30, 30)mm,  respectively.  Gaussian  noise  of  5%  was  added  to  the  simu¬ 
lated  light  intensity  change  on  the  detector  plane.  OPTICA  operates  on  a  noisy  scattering  wave 
— 0aca(rd,  r,)[l  +  n(rrj,  r,)]  where  n(rd.  rs)  is  a  Gaussian  random  variable  of  an  standard  deviation 
0.05. 

Independent  component  analysis  of  the  perturbations  in  the  spatial  intensity  distributions  pro¬ 
vided  corresponding  independent  intensity  distributions  on  the  source  and  detector  planes.  ICA 
generated  independent  intensity  distributions  on  the  source  and  detector  planes  are  shown  in  Fig. 
2,  for  the  two  absorptive  inhomogeneities.  Locations  of  the  absorptive  objects  are  obtained  from  fit¬ 
ting  these  independent  component  intensity  distributions  to  those  of  the  diffusion  approximation  in 
a  slab  Eq.  (4)  by  the  least  square  procedure  Eq.  (7).  The  first  object  is  found  at  (50.0, 60.0, 20.0)mm 
and  the  second  one  at  (30.0, 30.0, 30.1)mm.  The  coordinates  of  both  objects  agree  to  within  0.1mm 
of  their  known  locations.  The  strengths  of  the  two  objects  are  qi  =  1.00  and  #2  =  0.99  respectively, 
with  an  error  not  greater  than  1%  of  the  true  values. 

[Figure  2  about  here.] 

3.2  Discrimination  between  absorptive  and  scattering  inhomogeneities 

In  the  second  example,  one  absorptive  object  of  absorption  strength  of  0.01  is  placed  at  (50, 60, 20)mm 
and  one  scattering  object  of  diffusion  strength  of  negative  unity  (corresponding  to  an  increase  in 
scattering  for  the  inhomogeneity)  is  placed  at  (30, 30, 30)mm  respectively.  5%  Gaussian  noise  was 
added  to  the  simulated  light  intensity  change  on  the  detector  plane. 

Fig.  3  shows  the  ICA  generated  independent  intensity  distributions  on  the  source  and  detector 
planes  and  the  least  square  fitting.  The  first  column  corresponds  to  the  absorptive  inhomogeneity. 
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The  second  through  fourth  columns  correspond  to  the  scattering  object  which  produces  one  pair 
of  centrosymmetric  and  two  pairs  of  dumbbell-shaped  virtual  sources  and  mixing  vectors.  The 
absorptive  inhomogeneity  is  found  to  be  at  (50.2, 60.3, 20. 2)mm  with  a  strength  q\  =  0.0101.  The 
scattering  object  produces  three  pairs  (one  centrosymmetric  and  two  dumbbell-shaped)  of  virtual 
sources  and  mixing  vectors  centering  around  the  position  (x,y)  =  (30, 30)mm  [see  the  second 
through  fourth  columns  in  Fig.  3].  The  dumbbell-shaped  virtual  source  or  mixing  vector  comprises 
one  bright  part  and  its  antisymmetric  dark  counterpart.  The  resolved  position  and  strength  of  the 
scattering  object  are  found  to  be  (30.0, 30.0, 30.0)mm  and  q2  =  —0.99,  (32.1, 32.4, 30. 2)mm  and 
q2  =  —0.96,  and  (31.3, 30.2, 27.1)mm  and  <72  =  —1-05,  respectively,  through  fitting  to  the  individual 
pair.  For  the  scattering  object,  the  best  result  is  obtained  from  the  fitting  to  the  first  pair  of 
centrosymmetric  virtual  source  and  mixing  vector  from  the  scattering  object.  Taking  the  position 
and  strength  of  the  scattering  object  to  be  that  from  fitting  the  centrosymmetric  virtual  source  and 
mixing  vector,  the  error  of  the  resolved  positions  of  both  objects  is  within  0.3mm  of  their  known 
locations.  The  error  of  the  resolved  strengths  of  both  objects  is  approximately  1%  of  the  true  values. 

[Figure  3  about  here.] 

3.3  Co-located  absorptive  and  scattering  inhomogeneities 

For  one  complex  inhomogeneity  which  is  both  absorptive  and  scattering,  the  two  pairs  of  dumbbell¬ 
shaped  virtual  sources  and  mixing  vectors  produced  by  its  scattering  abnormality  can  be  used  to 
obtain  its  scattering  strength.  By  subtracting  the  scattering  contribution  off  the  measured  scattered 
wave,  our  procedure  can  be  applied  again  to  the  cleaned  data  and  proceed  to  obtain  its  absorption 
strength.  The  third  example  considers  a  complex  inhomogeneity  at  (30, 30, 20)mm  with  strengths 
of  absorption  qi  =  0.01  and  diffusion  q2  =  1  (corresponding  to  a  decrease  in  scattering)  respectively. 
5%  Gaussian  noise  was  added  to  the  simulated  light  intensity  change  on  the  detector  plane. 

Fig.  4  shows  the  ICA  generated  independent  intensity  distributions  on  the  source  and  detec¬ 
tor  planes  and  the  least  square  fitting.  The  first,  and  second  columns  correspond  to  the  pairs  of 
dumbbell-shaped  virtual  sources  and  mixing  vectors  produced  by  its  scattering  component.  The  po¬ 
sition  and  strength  of  this  diffusive  component  is  obtained  to  be  (32.7, 33.0,  20.5)mm  and  q2  =  0.95, 
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and  (31.7, 30.1, 20.4)mm  and  q2  =  0.96  by  fitting  the  two  individual  dumbbell-shaped  pair  respec¬ 
tively.  The  position  and  strength  of  the  diffusive  component  is  found  to  be  (30.9, 30.9, 20.4)mm 
and  q2  =  0.95  if  both  dumbbell-shaped  virtual  sources  and  mixing  vectors  are  used  in  fitting.  The 
third  column  of  Fig.  4  corresponds  to  its  absorptive  component  obtained  by  first  removing  the  scat¬ 
tering  contribution  from  the  measured  scattered  wave.  The  depth  and  strength  of  the  absorption 
component  is  found  to  be  (30.8, 30.7, 32. 7)mm  and  q\  =  0.0091. 

The  error  in  positioning  the  scattering  component  is  less  than  1mm  and  the  error  of  the  resolved 
strength  of  the  scattering  strength  is  ~  5%.  The  errors  in  positioning  and  the  resolved  strength  of 
the  absorptive  component  equal  to  ~  3mm  and  ~  10%,  respectively,  which  are  larger  than  those  of 
the  scattering  component  because  the  error  is  amplified  when  the  estimated  scattering  component 
is  used  in  subtraction  off  its  contribution  to  the  scattered  wave  in  our  procedure. 

[Figure  4  about  here.] 


3.4  Effect  of  noise 

To  demonstrate  the  effect  of  noise  on  the  performance  of  OPTICA,  different  levels  of  Gaussian  noise 
were  added  to  the  simulated  light  intensity  change  on  the  detector  plane. 

Fig.  5  shows  the  case  presented  in  Fig.  3  of  Sec.  3.2  with  instead  10%  and  20%  Gaussian  noise 
added  to  the  scattered  wave.  The  resolved  absorptive  inhomogeneity  is  at  (50.2, 60.3, 20.1)mm  with 
strength  0.0101  at  10%  noise,  and  at  (50.1, 60.3,  20.5)mm  with  strength  0.0102  at  20%  noise.  The 
resolved  position  and  strength  of  the  scattering  object  are  found  to  be  (30.0,  30.1, 30.0)mm  and  q2  = 
—0.98,  (32.1, 32.4, 30.4)mm  and  q2  =  —0.95,  and  (31.4, 30.1, 27.5)mm  and  <72  =  —1-00  respectively 
through  fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of  dumbbell-shaped  virtual  sources 
and  mixing  vectors  [see  the  second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  10%  noise.  The 
resolved  values  become  (28.9, 27.0, 32.9)mm  and  q2  =  —0.59  from  fitting  the  pair  of  centrosymmetric 
virtual  source  and  mixing  vector  [see  the  second  column  of  Fig.  5(b)],  and  (30.3, 32.3, 26.6)mm  and 
q2  =  —1.33  from  fitting  the  first  pair  of  dumbbell-shaped  virtual  source  and  mixing  vector  [see  the 
third  column  of  Fig.  5(b)],  respectively,  at  20%  noise.  The  dumbbell-shaped  virtual  source  on  the 
source  plane,  of  the  second  pair  of  dumbbell-shaped  virtual  source  and  mixing  vector,  is  deformed 
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and  the  fitting  is  not  shown  [see  the  fourth  column  of  Fig.  5(b)].  The  deformation  of  dumb  bell 
appears  first  on  the  source  plane  with  the  increase  of  noise  as  the  grid  spacing  on  the  source  plane 
is  larger  than  that  on  the  detector  plane  in  the  simulation. 

The  error  in  localization  and  characterization  of  scattering  inhomogeneities  increases  rapidly 
with  the  increase  of  noise,  from  ~  0.1mm  in  positioning  and  ~  2%  in  strength  at  10%  noise  to 
~  3mm  in  positioning  and  ~  50%  in  strength  at  20%  noise.  On  the  other  hand,  the  effect  of  noise 
on  localization  and  characterization  of  absorptive  inhomogeneities  is  much  smaller,  the  errors  at 
both  noise  levels  are  less  than  0.5mm  in  positioning  and  ~  2%  in  strength.  The  results  in  Sec.  3.2 
and  this  section  are  summarized  in  Tab.  1. 

[Figure  5  about  here.] 

[Table  1  about  here.] 

3.5  Effect  of  uncertainty  in  background 

In  the  examples  discussed  above,  we  have  assumed  the  light  intensities  change  measured  on  the 
detector  plane  is  obtained  with  an  exact  knowledge  about  the  background.  To  examine  the  effect  of 
uncertainty  in  background  optical  property  on  the  performance  of  OPTICA,  we  model  the  error  in 
the  estimation  of  the  background  absorption  or  diffusion  coefficients  as  a  uniform  Gaussian  random 
field  /( r).  The  Gaussian  noise  addressed  in  Sec.  3.4  is  set  to  be  zero  here.  OPTICA  operates  on  a 
“dirty”  scattering  wave  —  </>Sca(rd,  rs)  +  <5(/>sca(r<i,  rs)  where  5^sca(r<i,  rs)  is  the  change  in  the  scattered 
wave  from  that  of  a  uniform  background  of  absorption  fia  (or  diffusion  D)  to  that  of  a  background 
of  absorption  fia  +  /( r)  (or  diffusion  D  +  /( r)).  The  magnitude  of  the  background  uncertainty  is 
represented  by  the  signal  to  noise  ratio  (SNR)  defined  by 

OMD,,D,  lnl  Er(j£rJ^  sea  (rd>rs)|2 

SNR(dB)  =  101og10  - -2-  (13) 

Erd  Ers  l^sca^di  r.s) 

Figs.  6(a-c)  show  the  case  presented  in  Fig.  3  of  Sec.  3.2  with  40dB,  34dB  and  lOdB  SNR  due 
to  background  absorption  uncertainty,  respectively.  The  resolved  absorptive  inhomogeneity  is  at 
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(50.2, 60.3, 20. l)mm  with  strength  0.0101  at  40dB  SNR,  and  at  (50.2, 60.3, 20. l)mm  with  strength 
0.0100  at  34dB  SNR,  and  at  (50.6, 60.3, 20.3)mm  with  strength  0.0090  at  lOdB  SNR. 

The  resolved  position  and  strength  of  the  scattering  object  are  found  to  be  (30.1, 30.1, 30. 0)mm 
and  q2  —  —0.99,  (32.1, 32.9, 30.0)mm  and  q2  =  —0.95,  and  (31.4, 30.0, 27.5)mm  and  q2  =  —1.01 
respectively  through  fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of  dumbbell-shaped  virtual 
sources  and  mixing  vectors  [see  the  second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  40dB 
SNR.  The  resolved  values  become  (31.6, 31.7, 25.3)mm  and  q2  —  —0.52  and  (31.7,29.6,31.7)  and 
q2  =  —0.78  at  34dB  and  lOdB  SNRs,  respectively,  from  fitting  the  pair  of  centrosymmetric  virtual 
source  and  mixing  vector  [see  the  second  columns  of  Figs.  5(b), (c)].  The  dumbbell-shaped  virtual 
sources  and  mixing  vectors,  esp  dumbbell-shaped  virtual  sources  on  the  source  plane,  are  deformed 
and  the  fitting  are  not  shown  [see  the  third  and  fourth  columns  of  Fig.  5(b), (c)].  The  results  for 
the  influence  of  background  absorption  uncertainty  on  the  performance  are  summarized  in  Tab.  2. 

[Figure  6  about  here.] 

[Table  2  about  here.] 

Figs.  7(a,b)  show  the  case  presented  in  Fig.  3  of  Sec.  3.2  with  34dB  and  lOclB  SNR  due  to 
background  scattering  uncertainty.  The  resolved  absorptive  inhomogeneity  is  at  (50.1, 60.3, 20.1)mm 
with  strength  0.0100  at  34dB  SNR  and  at  (49.9,60.5, 20.1)mm  with  strength  0.0099  at  lOdB  SNR. 

The  resolved  position  and  strength  of  the  scattering  object  are  found  to  be  (30.0, 30.1, 30. 0)mm 
and  q2  —  —0.99,  (32.2, 33.0,  30. 0)mm  and  q2  =  —0.96,  and  (32.3,  29.3,  27.1)mm  and  q2  —  —1.08 
respectively  through  fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of  dumbbell-shaped  virtual 
sources  and  mixing  vectors  [see  the  second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  34dB 
SNR.  The  resolved  position  and  strength  of  the  scattering  object  are  found  to  be  (31.7, 31.1,  32.5)mm 
and  q2  =  —0.75,  (30.9,  31.4,  27.5)mm  and  q2  =  —1.08,  respectively  through  fitting  to  the  pair  of 
centrosymmetric  and  the  first  pair  of  dumbbell-shaped  virtual  sources  and  mixing  vectors  [see  the 
second  to  fourth  columns  of  Fig.  5(b)]  at  lOdB  SNR.  The  results  for  the  influence  of  background 
scattering  uncertainty  on  the  performance  are  summarized  in  Tab.  3. 

[Figure  7  about  here.] 
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[Table  3  about  here.] 


The  uncertainty  in  the  background  absorption  or  diffusion  coefficients  affects  the  performance 
of  OPTICA  in  a  similar  fashion  as  the  noise  does  discussed  in  Sec.  3.4.  The  error  in  localization 
and  characterization  of  scattering  inhomogeneities  increases  rapidly  while  the  error  in  localization 
and  characterization  of  absorptive  inhomogeneities  only  increases  mildly  with  the  increase  of  the 
uncertainty  in  the  background  optical  property.  The  uncertainty  in  background  scattering  has  a 
less  adverse  effect  on  the  performance  of  OPTICA  than  that  in  background  absorption. 

4  Discussion 

The  simulational  study  of  OPTICA  presented  in  this  article  demonstrates  its  potential  in  optical 
imaging  of  objects  in  turbid  media.  It  is  shown  to  be  able  to  locate  and  characterize  absorptive 
and  scattering  inhomogeneities  within  highly  scattering  medium.  In  particular,  OPTICA  can  dis¬ 
criminate  between  absorptive  and  scattering  inhomogeneities  and  locate  and  characterize  complex 
inhomogeneities  which  is  both  absorptive  and  scattering.  The  accuracy  of  localization  and  char¬ 
acterization  of  inhomogeneities  is  high.  In  the  cases  investigated  for  concentrated  inhomogeneities 
within  a  tissue  emulating  slab  of  thickness  of  50mm,  the  errors  in  resolved  object  locations  are  not 
greater  than  1mm  and  the  errors  in  the  resolved  optical  strengths  are  ~  2%  under  favorable  noise 
levels  and  reliable  background  estimations. 

Noise  at  higher  levels  and/or  larger  uncertainty  in  the  optical  property  of  the  background 
medium  makes  it  difficult  to  discriminate  between  absorptive  and  scattering  inhomogeneities.  In 
such  a  situation,  other  corroborative  evidences,  such  as,  multi-wavelength  measurements  are  re¬ 
quired  to  determine  the  nature  of  inhomogeneities.  Noise  at  higher  levels  and/or  larger  uncertainty 
in  the  optical  property  of  the  background  medium  also  introduces  larger  errors  in  localization  and 
characterization  of  scattering  inhomogeneities.  The  accuracy  of  localization  and  characterization 
of  absorptive  inhomogeneities  is  only  affected  mildly  by  the  amount  of  noise  and/or  uncertainty  in 
the  range  investigated. 

OPTICA  has  several  salient  features.  First,  OPTICA  provides  the  independent  components 

due  to  the  inhomogeneities  with  minimal  processing  of  the  data  and  does  not  have  to  resort  to  any 
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specific  light  propagation  model  for  obtaining  this  information.  Specific  light  propagation  models 
are  needed  only  in  the  later  stage  to  determine  location  and  optical  strength.  Second,  different 
geometries,  or  even  an  arbitrary  shaped  boundary,  can  be  used  with  OPTICA.  Although  we  used 
the  slab  geometry  in  the  work  reported  in  this  article,  the  approach  does  not  depend  on  any  specific 
geometry.  Third,  the  approach  is  fast  and  is  amenable  to  near  real  time  detection  and  localization 
of  objects  in  a  turbid  medium,  which  is  a  key  consideration  for  in  vivo  medical  imaging. 

As  it  is  well  known,  the  diffusion  approximation  to  RTE  which  is  widely  used  in  inverse  image 
reconstruction,  does  not  apply  when  the  separation  between  any  two  of  the  source,  the  inhomogene¬ 
ity  and  the  detector  is  small,  or  when  there  are  clear  regions  in  the  medium.  A  special  treatment  is 
also  required  when  the  medium  has  aligned  microstructures,  such  as,  myofibrils,  axons,  or  collagen 
fibers  in  tissues.  [43]  The  fact  that  a  prior  assumption  of  a  specific  model  of  light  propagation  in  the 
medium  is  not  assumed  in  the  identification  of  independent  components  by  ICA  and  only  required 
in  a  Green’s  function  analysis  of  the  retrieved  independent  component  is  desirable,  esp.  in  such 
situations  which  demands  a  more  complex  model  than  the  conventional  DA.  Performing  the  fitting 
procedure  for  each  identified  independent  component  is  much  simpler  and  transparent  than  match¬ 
ing  the  measured  light  intensity  to  a  forward  model  iteratively.  The  quality  of  reconstruction  of 
OPTICA  is  expected  to  be  higher  than  the  conventional  approach  when  only  an  imperfect  forward 
model  is  available. 

OPTICA  is  most  suited  to  detect  small  objects.  Given  its  ability  to  identify  low-contrast  small 
objects,  the  approach  is  expected  to  be  especially  useful  for  detection  of  tumors  at  their  early  stages 
of  development. 
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Multiple  passages  of  light  through  an  absorption  inhomogeneity  of  finite  size  deep  within  a  turbid  medium 
are  analyzed  for  optical  imaging  by  use  of  the  self-energy  diagram.  The  nonlinear  correction  becomes  more 
important  for  an  inhomogeneity  of  a  larger  size  and  with  greater  contrast  in  absorption  with  respect  to  the  host 
background.  The  nonlinear  correction  factor  agrees  well  with  that  from  Monte  Carlo  simulations  for  cw  light. 
The  correction  is  approximately  50% -75%  in  the  near  infrared  for  an  absorption  inhomogeneity  with  the 
typical  optical  properties  found  in  tissues  and  five  times  the  size  of  the  transport  mean  free  path.  ©  2004 
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The  main  objective  of  optical  imaging  of  turbid  media 
is  to  locate  and  identify  the  embedded  inhomogeneities 
by  essentially  inverting  the  difference  in  photon 
transmittance  in  the  time  or  frequency  domains  due 
to  the  presence  of  these  inhomogeneities.1-4  The  key 
quantity  involved  is  the  Jacobian,  which  quantifies 
the  influence  on  the  detected  signal  due  to  the  change 
of  the  optical  parameters  of  the  medium.  The  linear 
perturbation  approach  is  suitable  for  calculating  the 
Jacobian  for  only  a  small  and  weak  absorption  in¬ 
homogeneity  and  is  not  valid  when  the  absorption 
strength  is  large.6  This  failure  can  be  attributed  to 
the  multiple  passages  through  the  abnormal  site  by 
the  photon.  The  most  important  correction  is  the 
self-energy  correction,6  which  takes  into  account  the 
repeated  visits  made  by  a  photon  through  the  site  up 
to  an  infinite  number  of  times.  The  presence  of  other 
inhomogeneity  islands  can  be  ignored  because  the 
photon  propagator  decreases  rapidly  with  the  distance 
between  two  separate  sites. 

In  this  Letter  the  nonlinear  correction  for  an  absorp¬ 
tion  inhomogeneity  of  a  large  strength  due  to  repeated 
visits  by  the  photon  is  modeled  by  a  nonlinear  correc¬ 
tion  factor  (NCF)  to  the  linear  perturbation  approach. 
The  NCF  as  a  function  of  the  size  and  the  strength 
of  the  inhomogeneity  is  estimated  by  use  of  the  self- 
energy  diagram.  The  NCF  is  obtained  from  the  cu- 
mulant  approximation  to  the  radiative  transfer  and 
verified  by  Monte  Carlo  simulations  for  cw  light.  The 
magnitude  of  the  NCF  is  0.5-1  for  an  absorptive  in¬ 
homogeneity  of  up  to  5 lt  {It  is  the  mean  transport  free 
path  of  light)  and  of  the  typical  optical  properties  of 
human  tissues  (/ ualt/c  ~  0.01-0.05,  where  pta  is  the 
absorption  coefficient  and  c  is  the  speed  of  light  in  the 
medium). 

If  we  consider  an  absorption  site  centered  at  r  and 
far  away  from  both  the  source  and  the  detector,  change 
in  detected  light  A/  at  detector  rd  from  a  modulated 
point  source  at  rs  including  the  multiple  passages 
through  the  site  is  given  by 

0146-9592/04/151757-03$15.00/0 


A7  =  -G(rd,  to  IrjFS^Cr)  £  [-NscU(<o;R)V8pta(r)T 

0 

X  G(r,  u>  |rs) 

=  -G(rd,  (x>  |  r)  — — - — — — - — — 

1  +  Nse]f(to;  R)V8/i.a(r) 

X  G(r,  w|rs),  (1) 

where  8  /xa  is  the  excess  absorption  of  the  absorption 
site  of  size  R  and  volume  V,  to  is  the  modulation  fre¬ 
quency  of  light,  G  is  the  propagator  of  photon  migra¬ 
tion  in  the  background  medium,  and 

AW  («;#)=  J  Jv  Gfa,  to  |  ri)d’!r2d3ri  (2) 

is  the  self-propagator  that  describes  the  probability 
that  a  photon  revisits  volume  V.  Here  G(r2,  m  |  ri) 
gives  the  probability  density  that  a  photon  leaves  the 
volume  at  iq  and  reenters  it  at  r2.  The  scattering 
property  of  the  site  is  the  same  as  that  of  the  back¬ 
ground.  InEq.  (1)  G( rd,  co  |  rj  and  G(r,  to  |  rs)  are  well 
modeled  by  the  center-moved  diffusion  model  as  long 
as  separations  \rd  —  r|  and  |r,  —  r|  are  much  greater 
than  If.1  However,  the  diffusion  Green’s  function  can¬ 
not  be  used  in  Eq.  (2)  to  evaluate  Nse\[{co;R)  because 
the  diffusion  approximation  breaks  down  when  iq  is  in 
the  proximity  of  r2. 

Comparing  Eq.  (1)  with  the  standard  linear  pertur¬ 
bation  approach,  the  nonlinear  multiple  passage  effect 
of  an  absorption  site  is  represented  by  a  NCF: 

NCF  =  [1  +  Fseiff^fljVS/Ur)]-1.  (3) 

This  factor  serves  as  a  universal  measure  of  the  non¬ 
linear  multiple-passage  effect  as  long  as  the  absorption 
site  is  far  from  both  the  source  and  the  detector  and  its 
size  is  much  smaller  than  its  distance  to  both  the  source 
and  the  detector.  This  correction  is  more  significant 
when  the  NCF  is  further  away  from  unity. 

©  2004  Optical  Society  of  America 


1758  OPTICS  LETTERS  /  Vol.  29,  No.  15  /  August  1,  2004 


Photon  propagator  7V(r2, 1 1  ri,  s),  the  probability 
that  a  photon  propagates  from  position  iq  with  propa¬ 
gation  direction  s  to  position  r2  in  time  t,  for  any  sepa¬ 
ration  between  ri  and  r2,  was  recently  derived7,8  in 
the  form  of  a  cumulant  approximation  to  the  radiative 
transfer. 

In  the  case  of  interest  in  which  the  absorption  site 
is  deep  inside  the  medium,  the  photon  distribution 
is  isotropic.  The  photon  propagator  is  simplified 
to  Ne{[(r,t)  =  AIeff(|r2  -  nl,  t),  which  is  obtained  by 
averaging  AC(r2,  t  |  ri,  s)  over  the  propagation  direction 
s  of  light  over  the  47r  solid  angle.  In  the  frequency 
domain  this  effective  propagator  is  approximately 
given  by 


Net dr,  <o) 


w^exp 

exp(—xlt)  . 


+ 


4nDrKlt 


sinh(/cr) ,  r  <  l,’  (4) 


exp(-/cr)  .  w  , , 


where  D  =  ltc/ 3  and  k  =  [3(/ia  —  io))/ltc]1/2,  whose 
sign  is  chosen  with  a  nonnegative  real  part.  The  two 
terms  in  ATeff  when  r  <lt  represent  ballistic  and  dif¬ 
fusion  contributions,  respectively.  The  ballistic  term 
does  not  depend  on  scattering  because  the  photon  dis¬ 
tribution  involved  is  already  isotropic.  Only  diffusion 
contributes  to  when  r  >  lt.  The  self-propagator 
for  an  absorption  sphere  deep  inside  the  medium  is 
given  by 


N self  (w ;  R) 


V2  fv  fv  Neff^r2  ~  ril>")d3r2d3ri 


J_  f2» 

V  Jo  ‘ 


iVeff (r,  w)y0(r)4irrzdr , 


(5) 


where  yo{r)  =  1  -  (3r/4R)  +  {1/16)  (r/R)3  is  the 
characteristic  function  for  a  uniform  sphere.9  An  ab¬ 
sorption  site  of  an  arbitrary  shape  can  be  treated  the 
same  way.  The  exact  self-propagator  must  be  com¬ 
puted  by  a_numerical  quadrature.  A  good  approxi¬ 
mation  of  N se\f  {(o;R)  is 


Nseii(u>]R)  = 


Vc 


X 


(|f  +  e)-i3xit  +  &{K2),  f^i/2 

(-tz  +  -  -  —  r1  +  —  r3)  ’  (6) 

\5  e  2  16 &  320  g  J 

-£3Klt  +  0(k2),  £  >  1/2 


by  use  of  relation  (4),  where  £  =  R/l,  when  |k-|I?  1. 

The  exact  and  approximate  versions  of  dimensionless 
self-propagator  A/rseifW<"1c  when  k  =  0  are  plotted  as 
solid  and  dashed  curves,  respectively,  in  Fig.  1(a).  Di¬ 
mensionless  self-propagator  Ntie\(Vl/1c  depends  solely 
on  two  dimensionless  quantities  Klt  of  the  background 
and  R/l/  of  the  absorbing  sphere. 

It  is  worthwhile  to  point  out  that  the  self-propagator 
in  time  NS0]f(t-,R),  the  inverse  Fourier  transform  of 


Eq.  (5),  includes  the  contribution  from  the  ballistic  mo¬ 
tion  of  the  photon  when  the  photon  passes  through  the 
site.  This  ballistic  contribution  manifests  itself  as  the 
linear  decay  of  Nse.if{t;R)V  in  the  form  of  yo(ct)  near 
the  origin  of  the  time,  followed  by  a  transition  to  dif¬ 
fusion  [Fig.  1(b)]. 

The  NCF  is  obtained  by  plugging  Eq.  (5)  or  (6)  into 
Eq.  (3).  In  particular,  we  have 


NCF  = 


£  —  1/2 


>  (7) 


I  >1/2 


where  q  =  V 8 fia(r)/lfc  is  the  dimensionless  strength 
of  the  absorber  when  |/c|.R  «  1.  For  an  absorber  of 
fixed  q  >  0,  the  effectiveness  of  absorbing  light  is  di¬ 
minished  (the  NCF  decreases)  when  its  size  is  reduced. 
This  can  be  understood  from  the  fact  that  the  photon 
spends  less  time  per  volume  inside  the  absorber  of  a 
smaller  dimension  because  of  the  ballistic  motion  of  the 
photon  after  each  scattering  event.  The  photon  leaves 
a  small  site  (R  <  lt)  in  an  almost  straight  line.  The 
diffusion  behavior  for  an  individual  photon  is  observed 
only  after  a  large  number  of  scattering  and  on  a  scale 
larger  than  lt . 

Figure  2  shows  plots  of  the  NCF  versus  absorber 
size  for  typical  absorbers  of  excess  absorption  8/xalt/c 
equal  to  0.01  and  0.05.  The  nonlinear  correction  fac¬ 
tor  generally  decreases  with  the  size  of  the  absorber 
whose  excess  absorption  is  fixed.  With  the  increase 
of  the  background  absorption  and  the  modulation  fre¬ 
quency,  the  nonlinear  correction  becomes  less  accentu¬ 
ated.  The  phase  delay  is  larger  for  higher  modulation 
frequencies  and  less  background  absorption. 

Monte  Carlo  simulations10  are  performed  for  cw  light 
propagating  in  a  uniform  nonabsorbing  and  isotropic 
scattering  slab.  The  thickness  of  the  slab  is  L  =  80/; . 
A  spherical  absorber  of  radius  R  is  located  at  the  center 
(0, 0,  L/ 2)  of  the  slab.  The  excess  absorption  of  the  ab¬ 
sorber  is  8/j.alt/c  =  0.01.  The  absorber  has  the  same 
scattering  property  as  the  background.  The  details  of 
the  Monte  Carlo  computation  were  provided  in  a  pre¬ 
vious  publication.11  The  correlated  sampling  method 


Fig.  1.  (a)  Self-propagator  Nscif((y;.R)VZ(~1c  and  its  ap¬ 
proximation  form  when  k  =  0.  (b)  Self-propagator  for 

spheres  of  various  radii  in  the  time  domain  inside  a 
nonabsorbing  medium. 
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Fig.  2.  NCF  (magnitude  and  phase  angle)  versus  the  size 
of  absorbers  whose  excess  absorption  Sfj.alt/c  equals  0.01 
and  0.05.  Note  that  k212  =  3(/xa  -  i<o)lt/c  for  the  back¬ 
ground  medium. 


Fig.  3.  (a)  Theoretical  nonlinear  correction  factors  from 

numerical  quadrature  (Exact),  the  approximate  form  of 
relation  (7)  (Approx),  and  Monte  Carlo  simulations  (MC). 
Results  from  four  independent  Monte  Carlo  simulations 
are  shown  for  each  radius.  The  standard  linear  perturba¬ 
tion  approach  corresponds  to  horizontal  line  NCF  =  1  (not 
shown  in  the  figure),  (b)  Percentage  change  of  the  cw 
transmittance  from  the  experimental  data  given  in  Fig.  9 
of  Ref.  5  compared  with  the  theoretical  predictions  made 
by  the  standard  linear  perturbation  approach  (StdPert) 
and  those  including  NCF  (Exact  and  Approx). 


is  used  in  each  simulation  to  reduce  variance.12  A 
single  simulation  is  used  to  compute  the  change  in  light 
transmittance  due  to  the  presence  of  the  absorption 
site  and  the  corresponding  NCF.  Figure  3(a)  shows 
the  NCFs  obtained  from  numerical  quadrature,  the  ap¬ 
proximate  form  of  relation  (7),  and  Monte  Carlo  simu¬ 
lations.  The  agreement  between  our  theoretical  NCF 
and  that  from  Monte  Carlo  simulations  is  excellent. 
The  slight  difference  between  them  at  large  radii  is 
accounted  for  by  the  fact  that  the  sphere  can  no  longer 


be  regarded  as  small  compared  with  the  dimensions  of 
the  slab. 

Figure  3(b)  shows  the  percentage  change  of  the  cw 
transmittance  estimated  from  the  experimental  data 
given  in  Fig.  9  of  Ref.  5.  The  relevant  parameters 
of  the  experiment  are  summarized  in  the  inset.  The 
theoretical  predictions  from  the  linear  perturbation  ap¬ 
proach  with  and  without  the  nonlinear  correction  are 
also  shown  in  Fig.  3(b),  assuming  a  collimated  point 
source  and  a  point  detector  in  a  confocal  setup.  Our 
theoretical  prediction  with  nonlinear  correction  pro¬ 
vides  a  significant  improvement  over  linear  pertur¬ 
bation  and  agrees  much  better  with  the  experimental 
result. 

The  typical  value  of  the  absorption  coefficient  of 
human  tissues  in  the  near  infrared  indicates  that 
Hali/c  ~  0.01-0. 05. 13,14  This  fact  should  put  our 
results  on  NCFs  in  this  range  (Figs.  2  and  3)  into 
perspective.  The  nonlinear  correction  becomes  more 
important  for  an  inhomogeneity  of  a  larger  size  and 
with  greater  contrast  in  absorption  with  respect  to 
the  background.  The  value  of  the  NCF  decreases 
from  ~0.75  to  —0.5  for  an  absorption  site  of  radius  51, 
with  excess  absorption  SfA.alt/c  increasing  from  0.01  to 
0.05.  The  standard  linear  perturbation  approach  in 
optical  imaging  should  be  augmented  to  include  this 
nonlinear  correction. 
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Light-extinction  measurements  in  the  wavelength  range  of  400  to  1000  nm  are  performed  in  situ  on  Bacil¬ 
lus  subtilis  endospores  during  heat-shock-induced  activation.  Simultaneous  information  on  particle  size  and 
refractive  indices  during  activation  is  calculated  from  the  transmission  spectra  by  use  of  the  Gaussian  ray 
approximation  of  anomalous  diffraction  theory.  During  activation  the  refractive  index  of  the  core  decreases 
from  1.51  to  1.39,  and  the  size  increases  from  0.38  to  0.6  fim.  ©  2005  Optical  Society  of  America 
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Light  scattering  has  long  been  investigated  as  a  tool 
for  identifying  bacteria  size  and  shape,1-4  and  quasi¬ 
elastic  light  scattering  has  been  used  to  study  en- 
dospore  (ES)  structure.6,6  The  Gaussian  ray  approxi¬ 
mation  (GRA)  of  anomalous  diffraction  theory  has 
been  used  to  determine  the  size  and  shape  of  differ¬ 
ent  species  of  bacteria  from  light  transmission.3,7 
Bacterial  cell  size  typically  varies  from  the  submi¬ 
crometer  level  to  several  micrometers;  thus  their 
scattering  properties  in  the  visible  are  strongly  de¬ 
pendent  on  wavelength,  size,  shape,  and  refractive 
index.  The  bacteria  genera  Bacillus  and  Clostridium 
can  differentiate  to  form  an  ES,  a  dormant  cell  type, 
in  response  to  hostile  environments  and  can  revert  to 
vegetative  cells  by  germination  when  conditions  be¬ 
come  more  receptive.  The  structural  changes  occur¬ 
ring  during  germination  are  important  to  biologists5 
and  in  bioagent  detection.  In  the  past,  fluorescence 
was  shown  to  be  an  effective  technique  to  detect  ESs 
in  situ.s  Combining  light  scattering,  extinction,  and 
fluorescence  might  offer  an  ideal  tool  to  determine 
the  type  of  agent  present. 

In  this  Letter  light  extinction  is  used  in  situ  to 
monitor  Bacillus  subtilis  ES  activation  and  to  simul¬ 
taneously  retrieve  the  refractive  index  and  size  of  the 
ES  during  this  process.  The  optical  techniques  and 
methods  described  allow  for  the  real-time  visualiza¬ 
tion  of  the  ES  activation  process  and  the  initial  stage 
of  bacterial  differentiation.  The  results  may  have  im¬ 
portant  significance  in  bioagent  identification  sys¬ 
tems  and  modeling  applications. 

ESs  are  significantly  denser  than  and  their  refrac¬ 
tion  index  is  greater  than  that  of  vegetative  cells.  ESs 
consist  of  a  high-density,  dehydrated  core  surrounded 
by  a  lower-density  spore  coat  composed  of  cross- 
linked  polypeptides.  Germination  begins  with  activa¬ 
tion,  during  which  the  coat  is  shed  and  the  core  hy¬ 
drates  with  a  subsequent  decrease  in  density  and 
refractive  index  and  an  increase  in  size. 

From  the  recent  statistical  interpretation  of 
anomalous  diffraction  theory9,10  light  extinction  by 
soft  particles  is  mainly  determined  by  the  mean  and 
mean-squared  size  of  the  particles  in  the  GRA.  The 
scattering  cross  section  can  be  approximated  by 


C,  =  77 t1 


4n27r2 

2  (m  -  l)2(/z2  -  a2) 


4n47r4 

— {m  -  l)4(/z4  +  3<x4  +  G/jPo2) 


3X4 


(1) 


keeping  only  the  leading  two  terms  in  the  GRA  [Eq. 
(10)  in  Ref.  9],  Here  X  is  the  wavelength,  n  is  the  re¬ 
fractive  index  of  the  background  media,  m  is  the  rela¬ 
tive  refractive  index  of  the  particle,  fi={l)  is  the  mean 
light  path  through  the  scatterer,  and  a  is  given  by 
ct2=(/2)-(/)2.  The  first  term  on  the  right-hand  side  of 
Eq.  (1)  is  the  scattering  cross  section  in  the  interme¬ 
diate  case  limit.11,12  The  second  term  is  a  correction 
introduced  by  the  GRA  when  the  condition 
(2tt77X)|to-1|  <  <1  is  not  met.  In  the  absence  of  ab¬ 
sorption,  optical  extinction  K  is  given  by 


K  =  CgNL  =  77T2 


4n27T2 

— T~ar'1 
x2 


16n4774/3r4\ 

- 2 - 2VL,  (2) 

X4  12  / 


where  N  is  the  concentration  and  L  is  the  path 
length.  For  uniform  spheres,  a  and  (3  are  given  by 

460 

a=*2(m-l)2,  (m-1)4.  (3) 

81 

For  soft  particles  with  a  size  smaller  than  or  compa¬ 
rable  with  X  the  first  term  in  Eq.  (2)  is  sufficient  to 
describe  K,  which  is  linear  in  X-2.  For  particles  with  a 
larger  size  or  refractive  index  the  second  term  in  Eq. 
(2)  must  be  included  to  accurately  describe  K. 

In  the  experiment  B.  subtilis  ESs  (strain  NCTC 
3610)  were  isolated  from  their  cellular  debris,13  heat 
activated,14  and  placed  in  a  germination  medium.15 
The  germination  medium  had  a  minimal  carbon 
source  and  limited  amino  acid  content  to  restrict  cell 
growth  and  prevent  reversal  of  activation. 

Transmission  was  measured  from  400  to  1000  nm, 
a  region  in  which  the  bacteria  have  little  absorption, 
and  losses  are  due  to  scattering.  The  concentration 
was  1.03  X  108  ESs/ml  in  a  1-em  quartz  cuvette.  The 
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transmission  spectra  were  measured  with  a  halogen 
light  source  (HL2000,  Ocean  Optics,  Dunedin, 
Florida)  and  compact  spectrometer  (HR2000,  Ocean 
Optics),  both  coupled  to  optical  fibers.  The  spectra 
were  acquired  at  1-min  intervals  for  the  first  3.5  h  af¬ 
ter  heat  shock  and  at  30-min  intervals  between  3.5 
and  20  h.  Each  spectra  was  integrated  for  5  s. 

A  least-squares  fit  to  a  function  of  the  form 

K=C2\-2  +  C4\-i  (4) 

was  applied  to  each  transmission  spectrum.  Figure  1 
shows  K  plotted  as  a  function  of  X.-2  before  heat  shock 
and  at  representative  times  after  heat  shock:  t= 1, 
30,  60,  120,  and  360  min.  The  respective  least- 
squares  fits  are  also  plotted  in  Fig.  1.  Coefficients  C2 
and  C4  are  plotted  at  all  time  intervals  in  Fig.  2.  As 
seen  in  Fig.  1,  the  spectrum  from  the  spores  immedi¬ 
ately  after  heat  shock  is  nearly  identical  to  that  of 
the  spores  before  heat  shock,  indicating  that  the  re¬ 
sponse  of  the  spores  to  heating  is  not  instantaneous. 
It  can  also  be  observed  in  Figs.  1  and  2  that  the  larg¬ 
est  reduction  in  scattering  occurs  in  the  first  2  h,  pri¬ 
marily  caused  by  a  dramatic  decrease  in  the  refrac¬ 
tive  index  as  the  ESs  shed  their  coat  and  hydrate  the 
core.  From  2  to  3  h  the  rate  of  reduction  in  refractive 
index  slowed  while  the  radius  continued  to  increase. 
The  scattering  cross  section  remained  relatively  con¬ 
stant  from  3  to  8  h,  indicating  that  the  refractive  in¬ 
dex  and  spore  size  were  not  changing  significantly. 
After  8  h,  K  decreased  slowly  as  the  ESs  decayed,  re¬ 
leasing  their  contents  into  solution. 

With  a  uniform  spherical  model  for  the  ESs,  which 
neglected  the  contribution  of  the  coat,  the  radius  and 
index  were  calculated  from  C2  and  C4  with  Eq.  (2) 
and  expressions  (3).  It  was  observed  that  the  radius 
(plotted  in  Fig.  3)  remained  relatively  constant  in  the 
first  30  min  after  heat  shock.  The  ES  radius  was  also 
calculated  with  only  C2  and  making  the  assumption 
that  the  refractive  index  was  equal  to  that  of  vegeta¬ 
tive  cells,  1.39  (Refs.  3  and  16;  also  plotted  in  Fig.  3). 
The  radius,  calculated  by  both  methods  is  in  agree¬ 
ment  with  t>3  h,  confirming  that  by  3  h  the  ESs 
have  lost  their  coat  and  the  index  is  close  to  that  of 
vegetative  cells.  For  £<3  h  the  radius  calculated  by 


Fig.  1.  Light  extinction  as  a  function  of  wavelength  plot¬ 
ted  for  spores  before  heat  shock  and  at  t=l  min,  30  min, 
1  h,  2  h,  and  6  h  after  heat  shock  (dashed  curves).  The 
least-squares  fit  to  K=C2^~2+C4\~4  is  shown  for  each  spec¬ 
trum  (solid  curves). 


Time  (hours) 

Fig.  2.  Least-squares  fit  coefficients  C2  and  C4  plotted  as  a 
function  of  time. 

use  of  a  fixed  index  of  1.39  is  substantially  larger 
than  the  radius  calculated  from  C2  and  C4  and  incor¬ 
rectly  predicts  that  the  radius  is  decreasing  over 
time.  This  lack  of  agreement  demonstrates  that  the 
ES  refractive  index  during  the  early  time  period  is 
greater  than  1.39.  The  ES  radius,  calculated  from  C2 
and  C4,  is  0.35  /im  for  the  unshocked  spores  and 
0.34  /im  at  1  min  after  heat  shock.  The  radius  gradu¬ 
ally  increases  to  a  maximum  of  —0.6  fim  after  6  h.  In 
a  nutrient-rich  medium  the  ESs  would  continue  to 
germinate  and  grow,  eventually  reaching  the  size  of 
vegetative  cells.  The  refraction  index,  calculated  from 
Eq.  (2)  and  expressions  (3)  and  plotted  in  Fig.  4,  is 
1.55  for  both  the  unshocked  spores  and  at  1  min  after 
heat  shock.  The  index  decreases  steadily  in  the  first 
3  h,  after  which  it  remains  constant  at  1.39,  showing 
that  by  this  time  the  core  has  hydrated  and  become 
closer  to  vegetative  cells  in  density  and  composition. 

A  better  model  for  the  spore  structure  that  takes 
into  account  the  spore  coat  is  two  concentric  nonab¬ 
sorbing  spheres  consisting  of  a  thin  outer  shell — the 
spore  coat — that  slowly  dissolves  and  an  inner 
sphere — the  spore  core — that  changes  in  index  and 
size.  The  ratio  of  spore  coat  thickness  to  cell  radius  is 
designated  e.  The  initial  spore  coat  thickness  is  taken 
to  be  70  nm,17  and,  for  an  approximate  initial  radius 
of  400  nm,  e  is  initially  set  to  0.175.  Since  the  calcu¬ 
lated  spore  radius  (see  Fig.  3)  was  unchanged  during 
the  first  30  min,  our  model  assumes  that  the  spore 
coat  also  remains  unchanged  during  this  time.  This 
model  assumes  that  from  £  =  0.5  to  £  =  3  h  the  coat 
thickness  uniformly  decreases  from  its  initial  value 
to  zero.  Microscopy  confirmed  that  the  spore  coat 
completely  dissolved  within  3  h  of  heat  shock. 

In  the  GRA  for  concentric  spheres  the  geometric 
ray  inside  the  particle  may  bisect  both  the  core  and 
the  coat,  with  relative  refractive  indices  ra*  and  m0, 
respectively.  The  modified  geometric  path,  Z'=Zj(mj 
-l)+l0(m0-l),  is  given  by 

l’  =  2 (m0  -  l)(r02  -  h2)m  +  2(m,  -  m0){r2  -  h2)m,  h 

<  r  i 

=2(m0-l)(r20-h2)1/2,  h>rh  (5) 

with  the  distribution  function  for  h  given  by  p(h) 
=  (2 hlr2)  and  normalized  to  unity.  Evaluating  (l')  and 
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(l2)  with  Eq.  (5)  and  keeping  only  the  first-order 
terms  in  e  allows  a  and  j3  to  be  approximated  by 

a  ~  2(m,  -  l)2  -  8(TOj  -  1)( mt  -  m0)e  +  ... , 

460  544 

P  «  —{mi  -  l)4  -  —{mi  -  m0)[{mi  -  l)3 
ol  Z  ( 

-  2 (m,-  -  m0f]s  +  . . . ,  (6) 

where  m0  is  taken  to  be  1.39.  The  radii  and  refractive 
indices  were  recalculated  with  this  model  and  are 
also  plotted  in  Figs.  3  and  4,  respectively.  The  recal¬ 
culated  spore  radius  is  0.37  /im  before  heat  shock. 
For  the  first  30  min  after  heat  shock  the  radius  re¬ 
mains  essentially  unchanged  at  0.38  ^m  and  then 
slowly  increases.  The  refractive  index  is  1.515  for  the 
unshocked  spores,  which  is  in  agreement  with  the 
value  of  1.52  reported  by  Tuminello  et  al.4  After  heat 
shock  the  index  drops  linearly  to  1.39  over  a  2-h  pe¬ 
riod. 

Confocal  microscopy  finds  an  ellipsoidal  shape  for 
the  ESs  before  heat  shock  with  a  long  axis  of  1.25  yum 
and  short  axes  of  0.75  pm,  equivalent  in  light  extinc¬ 
tion  to  a  sphere  of  radius  0.43  yam,  which  is  in  agree¬ 
ment  with  our  concentric  sphere  model. 

In  conclusion,  the  GRA  has  been  used  to  measure 
refractive-index  and  size  changes  in  ESs  during  heat- 
shock  activation  from  light-extinction  measurements. 
Within  2  h  of  the  start  of  activation  the  refractive  in¬ 
dex  decreased  from  1.51  to  1.39.  The  radius  began  to 
increase  after  the  first  half  hour  from  0.38  to  0.6  yam 
at  4  h.  This  work  demonstrates  that  light  extinction 
can  monitor  in  situ  both  the  refractive  index  and  the 
size  of  ESs,  may  potentially  be  part  of  an  optical  bio¬ 
agent  detection  scheme,  and  can  be  used  for  modeling 
applications. 


Fig.  3.  ES  radius  plotted  as  a  function  of  time  after  heat 
shock.  The  radius  was  calculated  by  three  methods:  (1) 
from  C2  and  C4  assuming  a  uniform  sphere  (i.e.,  without  a 
spore  coat),  (2)  a  uniform  sphere  with  n  =  1.39,  and  (3)  two 
concentric  spheres  with  the  outer  sphere  (spore  coat)  dis¬ 
solving  as  described  in  the  text. 


Fig.  4.  ES  refractive  index  calculated  for  uniform  spheres 
(no  coat)  and  two  concentric  spheres  with  the  outer  sphere 
(spore  coat)  dissolving  as  described  in  the  text. 


This  work  was  supported  in  part  by  a  NASA  Uni¬ 
versity  Research  Center;  the  New  York  State  Office 
of  Science,  Technology  and  Academic  Research; 
Northrop  Grumman  Corp;  the  Charles  E.  Culpepper 
Biomedical  Pilot  Initiative  Grant;  and  the  U.S. 
Department  of  the  Army  (grant  DAMD17-02-1-0516). 
A.  Katz’s  e-mail  address  is  akatz@ccny.cuny.edu. 

References 

1.  A.  L.  Koch  and  E.  Ehrenfeld,  Biochim.  Biophys.  Acta 
165,  262  (1968). 

2.  P.  J.  Wyatt,  Nature  221,  1257  (1969). 

3.  A.  Katz,  A.  Alimova,  M.  Xu,  E.  Rudolph,  M.  Shah,  H.  E. 
Savage,  R.  Rosen,  S.  A.  McCormick,  and  R.  R.  Alfano, 
IEEE  J.  Sel.  Top.  Quantum  Electron.  9,  277  (2003). 

4.  P.  S.  Tuminello,  E.  T.  Arakawa,  B.  N.  Khare,  J.  M. 
Wrobel,  M.  R.  Querry,  and  M.  E.  Milham,  Appl.  Opt. 
36,  2818  (1997). 

5.  S.  E.  Harding  and  P.  Johnson,  Biochem.  J.  220,  117 
(1984). 

6.  A.  D.  Molina-Garcia,  S.  E.  Harding,  L.  de  Pieri,  N.  Jan, 
and  W.  M.  Waites,  Biochem.  J.  263,  883  (1989). 

7.  A.  Katz,  A.  Alimova,  M.  Xu,  H.  E.  Savage,  M.  Shah,  R. 
B.  Rosen,  and  R.  R.  Alfano,  Proc.  SPIE  4965,  73  (2003). 

8.  A.  Alimova,  A.  Katz,  H.  E.  Savage,  M.  Shah,  G.  Minko, 
D.  V.  Will,  R.  B.  Rosen,  S.  A.  McCormick,  and  R.  R. 
Alfano,  Appl.  Opt.  42,  4080  (2003). 

9.  M.  Xu,  M.  Lax,  and  R.  R.  Alfano,  Opt.  Lett.  28,  179 
(2003). 

10.  M.  Xu,  Appl.  Opt.  42,  6710  (2003). 

11.  P.  Chylek  and  J.  Li,  Opt.  Commun.  117,  389  (1995). 

12.  H.  C.  van  de  Hulst,  Light  Scattering  by  Small  Particles 
(Dover,  New  York,  1981). 

13.  D.  Dubnau,  Methods  Enzymol.  21,  430  (1971). 

14.  S.  J.  Foster  and  K.  Johnstone,  in  Regulation  of 
Procariotic  Development.  Structural  and  Functional 
Analysis  of  Bacterial  Sporulation  and  Germination,  I. 
Smith,  R.  A.  Slepecky,  and  P.  Setlow,  eds.  (American 
Society  for  Microbiology,  Washington,  D.C.,  1989),  p. 
89. 

15.  B.  Setlow,  E.  Melly,  and  P.  Setlow,  J.  Bacteriol.  183, 
4894  (2001). 

16.  M.  Jonasz,  G.  Fournier,  and  D.  Stramski,  Appl.  Opt. 
36,  4214  (1997). 

17.  A.  Driks,  Microbiol.  Mol.  Biol.  Rev.  63,  1  (1999). 


APPENDIX  7 


Fractal  Mechanisms  of  light  scattering  in  biological  tissue  and  cells 
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We  use  fractal  continuous  random  media  to  model  visible  and  near  infrared  light 
scattering  by  biological  tissue  and  cell  suspensions.  The  power  slaw  of  the  reduced 
scattering  coefficient,  the  anisotropy  factor  of  scattering,  and  the  phase  function 
are  derived.  The  fractal  fluctuation  of  the  dielectric  permittivity,  dependent  on  the 
fractal  dimension  Df  and  the  cutoff  correlation  length  Zmax,  is  shown  to  determine 
the  characteristics  of  light  scattering.  Good  agreement  with  experimental  results  are 
reported.  Implications  on  spectroscopic  tissue  diagnosis  are  discussed. 
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The  interaction  of  light  with  tissue  and  cells  is  the  underlying  mechanism  for  optical  biomedical 
technology  used  in  optical  imaging  and  spectroscopy  for  detection  of  pathology  changes.  The 
optical  properties  of  tissue  are  determined  by  chromophores,  microstructures  and  local  refractive 
index  variations.  Microstructures  in  biological  tissue  range  from  organelles  0.2— 0.5 fj,m  or  smaller, 
mitochondria  1  —  4 pm  in  length  and  0.3  —  0.7 yum  in  diameter,  nuclei  3  —  lO/rm  in  diameter,  to 
mammalian  cells  10  —  30 f.im  in  diameter.  The  refractive  index  variation  is  about  0.04  —  0.10 
for  biological  tissue  with  a  background  refractive  index  n0  =  1.35.'  Recently,  the  nature  of  light 
scattering  in  biological  tissue  has  been  actively  studied.1-6  As  many  biological  tissues  have  fractal¬ 
like  organization  and  are  statistically  self-similar,2,7-9  a  discrete  particle  model  of  scattering  centers 
in  tissue1,4  was  proposed  to  model  light  scattering  by  tissue.  The  discrete  particle  model  assumes 
that  the  refractive-index  variations  caused  by  underlying  microscopic  structures  can  be  treated  as 
spherical  particles  with  sizes  distributed  according  to  a  powerlaw: 

77(a)  =  r)0a3~Df  (1) 

where  77(a)  is  the  volume  fraction  of  spherical  particles  of  radius  0  <  a  <  amax  with  a  maximum 
radius  amax,  r]0  is  a  constant,  and  Df  is  the  fractal  dimension.  On  a  microscopic  scale  the  con¬ 
stituents  of  tissue  have  no  clear  boundaries  and  merge  into  a  quasi-continuum  structure.  Discrete 
particles  may  not  be  appropriate  to  describe  tissue  inhomogeneities.  As  the  refractive  index  varia¬ 
tion  in  biological  tissue  is  weak,  tissue  is  better  modeled  as  continuous  random  media  where  light 
scattering  is  not  due  to  the  discontinuities  in  refractive  index  but  rather  weak  random  fluctuations 
of  the  dielectric  permittivity.5 

In  this  Letter,  we  propose  to  use  fractal  continuous  random  media  to  model  light  scattering 
by  biological  tissue  and  cells.  The  correlation  function  R(r)  of  the  random  fluctuations  of  the 
dielectric  permittivity  depends  on  the  fractal  dimension  Df  and  the  cutoff  correlation  length  Zmax. 
Analytical  expressions  are  derived  for  the  power  law  of  the  reduced  scattering  coefficient,  the 


anisotropy  factor  of  scattering,  and  the  phase  function.  By  examining  the  existing  experimental 
results,  the  fractal  fluctuation  of  the  dielectric  permittivity  is  shown  to  determine  visible  and  near 
infrared  light  scattering  by  biological  tissue  and  cell  suspensions.  The  connection  of  the  proposed 
model  to  the  discrete  particle  model  and  implications  on  spectroscopic  tissue  diagnosis  will  also 
be  discussed. 

Assuming  tissue  is  statistically  space  homogeneous  and  isotropic,  the  correlation  function  of 
the  dielectric  permittivity  can  be  written  as  R(r)  =  (Se( r')Se(r'  +  r))  where  <5e  is  the  fluctuation 
of  the  dielectric  permittivity  from  the  background  value.  R( r)  is  proportional  to  the  correlation 
function  of  the  fluctuation  of  the  refractive  index  Rn(r )  and  <5c(r)  =  2nl(m  —  1)  where  m  is 
the  relative  refractive  index  at  position  r  when  the  fluctuation  of  the  refractive  index  is  weak. 
Light  scattering  by  the  continuous  random  medium  is  determined  by  the  power  spectrum  R  of  the 
random  fluctuations  of  the  dielectric  permittivity.  The  amplitude  scattering  matrix10  of  the  weakly 
fluctuating  continuous  random  medium  at  the  scattering  angle  9  is  given  by 
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where  k  =  27rn0/A  is  the  wave  number  and  A  is  the  wavelength  of  light  in  vacuum.  A  simple 
exponential  correlation  function  was  considered  by  Moscoso  et.  al.5  for  modeling  tissue  light 
scattering.  It,  however,  did  not  describe  essential  features  of  light  scattering  by  tissue  such  as  the 
power  law  of  the  reduced  scattering  coefficient. 

The  correlation  function  of  the  random  fluctuations  of  the  dielectric  permittivity,  in  the  fractal 
continuous  medium  model,  is  assumed  to  be  an  average  of  exponential  functions  weighted  by  a 
power  law  distribution  Eq.  (1)  for  the  correlation  length  l: 


where  e2  =  An^(m  -  l)2,  lmax  is  the  cutoff  correlation  length,  and  En(z )  =  exp (-z/t)tn~2dt. 
The  correlation  function  at  the  origin  is  R( 0)  =  e2Tv  where  Tv  is  the  total  volume  fraction  of 
scattering  centers  ( Tv  ~  0.3  in  soft  tissue).  The  correlation  function  decays  exponentially  as 
R(r)  =  e2rjolmJ?f  (j^)  exp  [~j~)  at  a  large  separation  r  lmax-  The  power  spectrum  is 
given  by 
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where  2Fi  is  the  Hypergeometric  function. 

The  amplitude  scattering  function  S(6)  in  Eq.  (2)  can  now  be  written  as 

nr  0  fU  F2k6l3 

|SW|^  =  k°-R(  2fcsi“  2>  =  L  2,[1  +  2(1  -WPP*1**  (5) 

where  fi  =  cos  6.  From  Eq.  (5),  the  anisotropy  factor  (mean  cosine  of  the  scattering  angle)  can  be 
found: 
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The  unnormalized  phase  function  is  given  by: 
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The  reduced  scattering  coefficient,  defined  as  ns(l  —  g ),  is  given  by 
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(Df  +  1  )(Df  -  1  )(Df  -  3)  25-D/  sin  (^t r) 
is  a  constant  dependent  only  on  Df.  We  have  assumed  /c/max  »  1  in  our  derivations  of  Eqs.  (7) 
and  (8). 

The  values  of  g  in  the  fractal  continuous  medium  model  versus  the  cutoff  size  parameter  klmax 
for  various  Df  are  displayed  in  Fig.  1(a).  g  is  larger  with  the  increase  of  lmax  and  the  decrease  of 
Df.  Light  scattering  Eq.  (5)  can  be  regarded  to  be  a  weighted  sum  from  components  of  different 
correlation  lengths  where  the  peak  contribution  occurs  at  kl*  =  klmax  when  Df  <  4  and  kl* 
may  be  different  from  klmax  when  Df  >  4.  The  value  of  kl*  when  Df  >  4  is  displayed  in 
Fig.  1(b).  The  power  law  of  //'  Eq.  (8)  is  obtained  thanks  to  the  cutoff  at  large  correlation  lengths 
as  g  approaches  unity  for  components  of  increasing  correlation  lengths.  By  fitting  experimental 
wavelength  dependence  of  fi's  to  Eq.  (8),  one  can  determine  both  the  fractal  dimension  Df  and  the 
parameter  e2g0. 

Comparisons  of  the  fractal  continuous  medium  model  for  biological  tissue  and  cells  to  experi¬ 
mental  results  are  in  order.  We  first  fit  the  theoretical  power  spectrum  Eq.  (4)  to  the  power  spectrum 
of  index  variations  in  mouse  liver  tissue  reported  by  Schmitt  et.  al.2  [see  Fig.  2(a)].  The  fitting 
is  excellent  and  yields  Zmax  =  2.3gm  and  Df  =  4.0  for  mouse  liver  tissue.  A  single  exponential 
correlation  function5  will  not  fit.  We  do  not  know  any  results  on  the  wavelength  dependence  of  the 
reduced  scattering  coefficient  for  mouse  liver  tissue.  The  wavelength  dependence  of  the  reduced 
scattering  coefficient  of  rat  liver  tissue  was  reported  by  Parsa  et.  al.11  The  reduced  scattering  co¬ 
efficients  within  the  range  600  to  1400nm  are  displayed  in  Fig.  2(b)  and  fitted  well  to  a  power 


law  A-094.  The  fractal  dimension  of  rat  liver  tissue  is  hence  Df  =  3.94.  The  cutoff  correlation 
length  is  found  to  be  lmax  =  1.5/im  from  g  =  0.94  of  the  rat  liver  tissue  at  800nm.  Evidently,  both 
the  fractal  dimension  and  the  cutoff  correlation  length  extracted  from  the  power  spectrum  of  the 
refractive  index  variations  of  mouse  liver  tissue  agree  reasonably  well  with  those  extracted  from 
the  light  scattering  data  for  rat  liver  tissue.  This  gives  a  strong  support  to  the  fractal  continuous 
medium  model  for  light  scattering  by  tissue. 

Fig.  (3)  displays  the  phase  function  for  suspensions  of  rat  embryo  fibroblast  cells  (Ml)  and 
mitochondria,  respectively,  reported  by  Mourant  et.  al.3  and  the  fitting  to  the  theoretical  phase 
function  (7).  The  fractal  dimension  is  found  to  be  Df  =  3.86  and  4.58  for  Ml  cells  and  mito¬ 
chondria,  respectively.  The  value  of  the  fractal  dimension  for  Ml  cells  was  4  ±  0.07  using  the 
reported  wavelength  dependence  g!s  oc  A-1'0*0  07  over  the  wavelength  range  500  —  800nm.  The 
agreement  between  the  two  values  of  Df  for  Ml  cells  from  fitting  either  the  phase  function  or 
the  powerlaw  of  g!s  is  good.  The  g-factor  for  Ml  cells  was  reported  to  be  0.98  and  the  maxi¬ 
mum  correlation  length  can  be  estimated  to  be  lmax  ~  3.2//,m.  The  ^-factor  for  mitochondria 
can  be  computed  from  the  phase  function  to  be  0.81  and  the  maximum  correlation  length  is  then 
estimated  to  be  lmax  ~  0.6//m.  The  component  contributing  most  to  light  scattering  has  size  pa¬ 
rameter  kl*  =  klmax  ~  40  and  kl*  ~  1  in  Ml  cells  and  mitochondria,  respectively.  Ml  cells  have 
much  larger  scattering  centers  and  much  smaller  fractal  dimension  than  mitochondria.  The  larger 
scattering  centers  in  Ml  cells  are  due  to  the  nucleus. 

These  comparisons  show  that  the  fractal  continuous  medium  model  describes  well  light  scatter¬ 
ing  from  both  biological  tissue  and  cell  suspensions.  Cell  suspension  can  be  regarded  as  composi¬ 
tion  of  many  mini-continuous  media  of  random  fluctuating  dielectric  permittivity  due  to  individual 
cells.  The  powerlaw  of  the  reduced  scattering  coefficient  originates  from  the  underlying  fractal 


fluctuation  of  the  refractive  index  of  the  medium. 


Light  scattering  in  the  fractal  continuous  medium  model  is  caused  by  weak  random  fluctuations 
of  the  dielectric  permittivity.  This  model,  however,  bears  a  close  connection  with  the  discrete 
particle  model.  We  can  approximate  the  amplitude  scattering  matrix  of  the  spherical  particles  in 
the  discrete  particle  model  by12 

2  =  2  |m  —  1 12  z6  (10) 

[1  +  2(1  —  cos0)x2] 

where  x  =  ka  is  the  size  parameter  of  the  particle  as  the  particles  are  soft  (|  m  —  1|  -C  1).  Ignoring 
the  effect  due  to  correlated  scattering  based  on  the  fact  that  correlated  scattering  among  particles 
of  size  much  less  than  a  wavelength  is  most  significant  while  their  contribution  to  scattering  is 
minimal13,  the  discrete  particle  model  assuming  a  particle  size  distribution  of  the  powerlaw  (1) 
reaches  the  same  amplitude  scattering  function  (5)  as  in  the  fractal  continuous  medium  model. 
This  illustrates  the  correlation  length  l  in  the  fractal  continuous  medium  model  may  be  intuitively 
linked  to  the  radius  of  “fictional”  scattering  centers  present  within  tissue. 

The  power  in  the  powerlaw  of  fi's  a  A-6  is  usually  called  the  scattering  power.  The  scattering 
power  in  the  fractal  continuous  medium  model  relates  to  the  fractal  dimension  of  the  underlying 
fluctuation  of  the  refractive  index  (b  =  Df  —  3)  and  should  be  distinguished  from  that  due  to  Mie 
particles  of  narrow  size  distribution.3  The  scattering  power  has  been  recognized  to  be  an  important 
parameter  in  discriminating  normal  and  cancerous  tissue.14-16  Both  the  fractal  dimension  Df  and 
the  parameter  e2rj0  can  be  estimated  from  fitting  the  wavelength  dependence  of  n's  to  Eq.  (8)  in  the 
fractal  continuous  medium  model.  The  value  of  Df  reveals  the  relative  weight  of  small  scattering 
centers  vs  large  scattering  centers.  The  value  of  e2r/()  represents  the  overall  density  of  scattering 
centers  which  is  proportional  to  the  radiographic  density  of  tissue  of  predictive  value  for  cancer 
risk.17  The  maximum  correlation  length  /max  can  be  estimated  from  the  anisotropy  factor.  The 
access  to  all  these  parameters  will  yield  much  more  valuable  information  about  the  structure  and 
the  physiological  state  of  tissue  than  using  the  scattering  power  alone. 


In  conclusion,  we  have  shown  that  light  scattering  properties  of  biological  tissue  and  cell  sus¬ 
pensions  can  be  well  represented  by  a  fractal  continuous  random  medium  model  where  light  scat¬ 
tering  is  due  to  weak  random  fluctuations  of  the  dielectric  permittivity.  The  fractal  dimension 
Df  and  the  cutoff  correlation  length  lmax  characterizes  the  essential  features  of  light  scattering  by 
biological  tissue  and  cell  suspensions,  including  the  wavelength  dependence  of  the  reduced  scat¬ 
tering  coefficient,  the  phase  function,  and  the  anisotropy  factor.  The  fractal  continuous  random 
medium  model  should  facilitate  the  analysis  of  light  scattering  spectroscopy  for  tissue  diagnosis 
and  provide  valuable  insight  to  light  scattering  mechanisms  by  tissue  and  cells. 

This  work  is  supported  in  part  by  NASA.  M.  X.  thanks  the  support  by  the  Department  of  Army 
(Grant#  DAMD 1 7-02- 1-0516). 
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Fig.  1  (a)  The  anisotropy  factor  of  the  fractal  continuous  random  medium  versus  the  cutoff  size 

parameter  klmax.  (b)  Size  parameter  kl*  of  the  component  contributing  most  to  light  scattering 
when  Df  >  4. 


Fig.  2  (a)  Power  spectrum  of  refractive  index  variations  in  mouse  liver  tissue  fitted  to  the 

theoretical  power  spectrum  Eq.  (4).  Symbols  represent  data  reported  by2  and  the  dash  line  is  the 
theoretical  fit.  Fitting  yields  Df  =  4.0  and  lmax  =  2.3//m.  (b)  The  wavelength  dependence  of 
the  reduced  scattering  coefficient  of  rat  liver  tissue  fitted  to  the  powerlaw  A 3~Df.  The  symbols 
represent  data  reported  by1 1  and  the  solid  line  show  the  fitted  curve.  Fitting  yields  Df  =  3.94  and 
^max  =  1.5/im. 


Fig.  3  Phase  function  of  suspensions  of  rat  embryo  fibroblast  cells  (M 1 )  and  mitochondria  fitted 
to  the  theoretical  phase  function  Eq.  (7).  Symbols  are  data  reported  by3  and  the  solid  lines  are 
theoretical  curves.  Fitting  yields  Df  =  3.86  and  4.58  for  Ml  cells  and  mitochondria  respectively. 
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